%% Interpolacja wielomianem
% opracował dr inż. Robert JAKUBOWSKI, Politechnika Rzeszowska, Katedra Samolotów i Silników Lotniczych
% Table of Contents
% Interpolacja wielomianem
% Wielomian interpolacyjny Lagrange'a
% Funkcje interpolacyjne w Matlabie
% Ćwiczenia do tematu
% Funkcje przygotowane do tematu
% Literatura
% Interpolacją nazywamy postępowanie prowadzące do znalezienia wartości pewnej funkcji f(x) w dowolnym punkcie przedziału (x0, xn) na podstawie znanych wartości tej funkcji w punktach x0, x1, ..., xn, zwanych węzłami interpolacji (x0 < x1 < ... < xn) [1].
% Postępowanie prowadzące do znalezienia wartości funkcji f(x) w punkcie leżącym poza przedziałem (x0, xn) nazywamy ekstrapolacją [1].
%
% Rys 1. Przykład inerpolacji i ex1trapolacji
% Najczęściej stosowaną metodą interpolacji jest interpolacja wielomianem
%% Interpolacja wielomianem
% Interpolacja wielomianem polega na znalezieniu wielomianu interpolacyjnego, który będzie przechodził dokładnie przez zadane punkty węzłowe. Zgodnie z teorią na bazie n punktów węzłowych jesteśmy w stanie opisać wielomian interpolacyjny co najwyżej n-1go stopnia. To znaczy na dwóch punktach jesteśmy w stanie opisać wielomian pierwszego stopnia, na trzech punktach drugiego stopnia na trzech - czwartego, na 50 punktach - 49-go stopnia (patrz rys. 2)
%
% Rys 2. Przykłady funkcji interpolacyjnej dla dwóch, trzech i czterech punktów.
% Zadaniu interpolacji sprowadza się do znalezienia wskazźników wielomianu funkcji interpolującej na bazie posiadanych punktów. Omówmy zagadnienie na podstawie n punktów.
% Dla zadanych n punktów czyli wektora x zawierającego n elementów i odpowiadającego temu wektora y zawierającego także n elemntów zadanie sprowadza się do wyznaczenie parametrów wielomianu n-1 stopnia czyli warości a(n-1) do a(0), czyli n-elementowy wektor a. Aby to zrobić potrzebne jest n równań. Te równania otrzymamy pisząc dla każedgo punktu równanie wiążace x1 poprzez zależność funkcyjną z wartością y dla danego pynktu.
% \begin{array}{l}
% a_n x_1^{n-1} +a_{n-1} x_1^{n-2} +\ldotp \ldotp \ldotp +a_3 x_1^2 +a_2 x_1 +a_1 =y_1 \\
% a_n x_2^{n-1} +a_{n-1} x_2^{n-2} +\ldotp \ldotp \ldotp +a_3 x_2^2 +a_2 x_2 +a_1 =y_2 \\
% \vdots \\
% a_n x_{n-1}^{n-1} +a_{n-1} x_{n-1}^{n-2} +\ldotp \ldotp \ldotp +a_3 x_{n-1}^2 +a_2 x_{n-1} +a_1 =y_{n-1} \\
% a_n x_n^{n-1} +a_{n-1} x_n^{n-2} +\ldotp \ldotp \ldotp +a_3 x_n^2 +a_2 x_n +a_1 =y_n
% \end{array}
% Sprowadza się to zatem do rozwiązania ukłądu równań liniowych ze względu na a_i, który w zapisie macieżowym można przedstawić:
% \left\lbrack \begin{array}{cccccc}
% x_1^{n-1}  & x_1^{n-2}  & \cdots  & x_1^2  & x_1  & 1\\
% x_2^{n-1}  & x_2^{n-2}  & \cdots  & x_2^2  & x_2  & 1\\
% \vdots  &  &  &  &  & \vdots \\
% x_{n-1}^{n-1}  & x_{n-1}^{n-2}  & \cdots  & x_{n-1}^2  & x_{n-1}  & 1\\
% x_n^{n-1}  & x_n^{n-2}  & \cdots  & x_n^2  & x_n  & 1
% \end{array}\right\rbrack *\left\lbrack \begin{array}{c}
% a_n \\
% a_{n-1} \\
% \vdots \\
% a_2 \\
% a_1
% \end{array}\right\rbrack =\left\lbrack \begin{array}{c}
% y_1 \\
% y_2 \\
% \vdots \\
% y_{n-1} \\
% y_n
% \end{array}\right\rbrack
% Rozwiązanie sprowadza się więc do przygotowania odpowiednich macierzy na bazie punktów x i y, a następnie jedną z poznanych metod rozwiązania układu równań liniowych. Przykład funkcji relizującej takie zadanie przedstawiono w rozdziale Funkcje przygotowane do tematu pod nazwą interpolacja
% Sprawdźmy działanie tej funkcji. Zadajmy pięć punktów w postaci wektora x1 i y1 i narysujmy je na wykresie:

x1=1:5;
y1=[0 2 2 1.5 5];
plot(x1,y1,"ro")

% Dla tak przedstawionych punktów wyznaczmy wartość wskaźników wielomianu interpolacyjnego z wykorzystaniem przygotowanej funkcji:

a=interpolacja(x1,y1)

% Otrzymane w ten sposób wartości wielomianu możemy przeliczyć dla dowolnego x na wartość y. Wykorzystamy istniejącą w Matlabie funkcję polyval. Funkcję polyval wywołuje się w postaci y=polyval(a,x), gdzie a jest wektorem określającym parametry wielomianu poczwszy od stojącego przy najwyższej potedze, do tego przy x^0, zaś x jest parametrenm dla którego ma być wyznaczona wartość y.Wyznaczmy zatem wartości funkcji dla określonych x1 i porónajmy z wartościami y1

y_1=polyval(a,x1)

% Otrzymane wartości y_1 odpowidają zadanym wartościom y1.
% Narysujmy funkcję w przedziale od x1(1) do x1(n) wykorzystując funkcję fplot

hold on
fplot(@(x) polyval(a,x),[x1(1),x1(end)],'g',"LineWidth",1.5)
grid on
hold off
xlabel("x")
ylabel("y")
legend ("punkty węzłowe", "fun  interpol");
hold off

%
%% Wielomian interpolacyjny Lagrange'a
% Najczęściej stosowanym podejściem do obliczeń wielomianu interpolującego jest metoda Lagrange'a. Wielomian interpolacyjny Lagranga otrzymuje się korzystając z zależności:
% W\left(x\right)=\sum_{j=1}^n \left(y_j \prod_{i=1\;i\;i\not= j}^n \frac{\left(x-x_i \right)}{x_j -x_i }\right)
% Przedstawmy zagadnienie dla punktów podobnie jak poprzednio

x1
y1

% Obliczenia wykonamy wykorzystując obliczenia symboliczne. Wprowadźmy zmienną symboliczną x1

syms x

% Rozpiszmy teraz posczególne składniki sumy wielomianu. Zgodnie z przedstawioną zależnością będzie ich tyle, co punktów wekściowych 5. Zgodnie z zależnością poszczególne składniki przedstawiją się:

W1(x)=y1(1)*(x-x1(2))*(x-x1(3))*(x-x1(4))*(x-x1(5))/...
    ((x1(1)-x1(2))*(x1(1)-x1(3))*(x1(1)-x1(4))*(x1(1)-x1(5)))
W2(x)=y1(2)*(x-x1(1))*(x-x1(3))*(x-x1(4))*(x-x1(5))/...
    ((x1(2)-x1(1))*(x1(2)-x1(3))*(x1(2)-x1(4))*(x1(2)-x1(5)))
W3(x)=y1(3)*(x-x1(1))*(x-x1(2))*(x-x1(4))*(x-x1(5))/...
    ((x1(3)-x1(1))*(x1(3)-x1(2))*(x1(3)-x1(4))*(x1(3)-x1(5)))
W4(x)=y1(4)*(x-x1(1))*(x-x1(2))*(x-x1(3))*(x-x1(5))/...
    ((x1(4)-x1(1))*(x1(4)-x1(2))*(x1(4)-x1(3))*(x1(4)-x1(5)))
W5(x)=y1(5)*(x-x1(1))*(x-x1(2))*(x-x1(3))*(x-x1(4))/...
    ((x1(5)-x1(1))*(x1(5)-x1(2))*(x1(5)-x1(3))*(x1(5)-x1(4)))

% Zsumujmy te skladniki w celu wyznaczenia wielomianu interpolacyjnego:

W(x)=W1+W2+W3+W4+W5

% Wykorzystując funkcję expand uprośćmy to wyrażenie:

W1(x)=expand(W)

% Sprwdźmy poprawność wyznaczonego wielomianu obliczając wartości dla przyjętych punktach węzłowych

y_2=W(x1);
y_2=double(y_2)

% Wyznaczone wartości y_2 są zgodne z wartościami y1, czyli w węzłach otrzymaliśmy takie same wartości, jak przyjęto do interpolacji, co wskazuje, że metoda działa poprawnie.
% Procedurę wyznacznia wielomianu interpolacyjnego metodą Lagrange'a zawarto w funkcji interpolacja_l. W wyniku działania funkcji otrzymuje się wynik będący funkcją zmiennej x. Sprwdżmy wynik działania funkcji dla wcześniej przyjętych punktów x1 i y1


syms x

WL(x)=interpolacja_l(x1,y)

%  Otrzymana funkcją jest dokładnie taka, jak poprzednio wyznaczono to na przykłądzie zrobionym wcześniej.
%
%% Funkcje interpolacyjne w Matlabie
%
%  Do interpolacji w Matlabie można wykorzystać funkcję polyfit. Na bazie podanych punktów pozwala ona wyznaczyć prametry wielomianu interpolacyjnego. Wywołanie funkcji P=polyfit(x,y,n). Zmiennymi wejściowymi są wektory x i y, które określają położenie punktów węzłowych oraz stopień wielomianu n. Dla interpolacji stopień n powinien mieć wartość n-1. Otrzymany wektor P zawira parametry wielomianu interpolacyjnego w postaci P=[P1, P2, P3 .. Pn], takie że
% W\left(x\right)={\mathrm{P}}_1 x^{n-1} +P_2 x^{n-2} +\ldotp \ldotp +P_{n-1} x+P_n
% Sprawdźmy działanie tej funkcji. Ponieważ daje ona takie same wyniki jak pierwsza z prezentowanych metod, gdzie wartości wskaźników wielomianu interpolacyjnego wyznaczyliśmy jako a. Za x podstawimy x1, a za y - y1. Ponieważ mamy wektor pięcioelementowy, dlatego stopień wielomianu interpolacyjnego jest n=
% Przypomnijmy wartości x1, y1 i a:

x1
y1
a

% Poszukajmy parametrów wielomianu interpolacyjnego metodą polyfit

P=polyfit(x1,y1,4)

% P jest takie same jak a, a zatem otrzymaliśmy dokładnie ten sam wielomian
% W tym przypadku podobnie jak poprzednio do wyznaczenia wartości funkcji interpolujacej dla danego x wykorzystamy funkcję polyval. Pokażmy działanie metody na przykładzie kilku punktów, które przedstawimy graficznie.

plot(x1,y1,'or')

% Wprowadzamy punkty dla których wyznaczymy wartości funkcji z wykorzystaniem wyznaczonego wielomianu

x2=x1(1)-1:0.3:x1(end)+0.5

% Wyznaczamy wartości y2

y2=polyval(P,x2)

% Wyznaczone wartości punktów przedstawimy graficznie punktami w kolorze zielonym

hold on
grid on
plot(x2,y2,'gp')
xlabel("x")
ylabel("y")
legend ("punkty węzłowe", "fun  interpol");
hold off

% Gdy nie interesuje nas wyznaczenie funkcji interpolującej, a jedynie wyznaczenie wartości pomiędzy węzłami interpolacyjnymi, to można skorystać z funkcji interp1. Wywołanie funkcji jest następujące y=interp1(x1,y1,x,metoda). Zmiennymi wejściowymi w tym przypadku jest wektor x dla którego chcemy wyznaczyć wartości y, oraz dane punktów węzłowych x1, y1. Można opcjonalnie wskazać metodę interpolacji metoda np "linear" "pchip" "spline" itp - inne metody patrz Matlab help.
% Przeanalizujmy jakie wyniki interpolacji dostaniemy wykorzystując tą funkcję z różnymi metodami. Weźmy pod uwagę te same węzły interpolacji, co wcześniej, a wartości x dla których wycznaczymy parametry wielomianu będą znajdowały się pomiędzy nimi. Na przykład przyjęte x_

x_=x1(1):0.2:x1(end);

% Wyznaczone wrtości y dla tej metody zapiszemy jako y_lin

y_lin=interp1(x1,y1,x_,"linear");
plot(x_,y_lin,'xg')
hold on
grid on

% Tym razem skorzystamy z interpolacji funkcją trzeciego stopnia  pomiędzy węzłami metoda cubic, obecnie rekomendowana pod nazwą pchip.

y_cub=interp1(x1,y1,x_,"pchip");
plot(x_,y_cub,"cs")

% Wartości wyznaczone metodą interpolacji splainem zapiszemy jako y_spl

y_spl=interp1(x1,y1,x_,"spline");
plot(x_,y_spl,"m^")
plot(x1,y1,'or','MarkerFaceColor',"r")
xlabel("x")
ylabel("y")
legend ("interpol lin", "interpol cub", "interpol spl", "punkty węzłowe","Location","northwest")
hold off

% W zależności od przyjętej metody obliczeń punktów pomiędzy węzłami w prezentowanej metodzie mamy rózne wartościdla punktów pośrednich
%% Ćwiczenia do tematu
% Utworzyć wektory x1 y1  zawierających osiem pónktów, Następnie dla tych punktów znaleźć wielomian interpolacyjny w postaci funkcji W(x)  wykorzystując funkcję interpolacja_l. W kolejnym kroku wyznaczyć wartości parametrów wielomianu interpolacyjnego inną metodą i porównać wyniki.
% Na koniec narysować punkty węzłowe, i wielomian na wykresie dodać legendę i opis osi.




%
%
%% Funkcje przygotowane do tematu
% Funkcja interpolacja do wyznaczania wskażników wielomianu interpolacyjnego dla zadanych punktów w postaci wektora x y. Zmienna wyjściowa a zawiera parametry wielomianu interpolacyjnego także a1 stoi przy najwyższej potędze x, a an jest wyrazem wolnym

function a=interpolacja(x,y)
% sprawdzenie ilości elementów wektora x i y
[s1,s2]=size(x);
[s3,s4]=size(y);
if or(s1~=s3,s2~=s4)
    % brak rozwiązania, ze względu na nierówną ilość elementów x i y
    a=NaN;
    display "Wymiar wektora x i y powinien być taki sam"
else
    
    if s1==1
        % odwracanie wektora w przypadku gdy jest on wektorem wierszowym
        y=y';
        x=x';
    end
    %================= W tej części jest zawarte rozwiązanie zadania
    if s1==1 || s2==1 %funkcja lub zpisana inaczej niż wczesniej
        n=length(x);
        i=0;
        %Tworzenia wektora x1 do rozwiązania ukąłdu równań x1*A=Y
        for j=n-1:-1:0
            i=i+1;
            x1(:,i)=x.^j;
        end    
        %Wyliczenie rozwiązania -rozwiązanie ukłądu równań ze względu na A
        a=x1^-1*y;
        % Zamiana a z kolumny na wiersz
        a=a';
    %==================================================
    else
        % brak rozwiązania, gdy zmienne są macierzami zamiast wektorów
        a=NaN;
        display "Obydwie zmienne powinny być wektorami o takiej samej ilości elementów"
    end
end    
end

% Funkcja interpolacja metodą Lagrange'a. Zmienne wejściowe, to zbiór współrzędnych węzłów interpolacyjnych x1 i y1. Wynikiem jest funkcja zmiennej x- wielomian interpolacyjny - W(x)

function W=interpolacja_l(x1,y1)
    n=length(x1);
    % wprowadzamy zmienną symboliczną x
    syms x
    %Przypisujemy początkową wartość wielomianu interpolacyjnego równą 0 
    W(x)=0*x;
    %Iteracyjnie w n krokach wyznacza się wielomian interpolacyjny
    for i=1:n
        % wprowadzamy wektor pomocniczy xp który zawiera wszystkie x1 za
        % wyjątkiem x1(i)
        xp=[x1(1:i-1),x1(i+1:n)];
        Wp(x)=1+0*x;
        
        for j=1:n-1
            Wp(x)=Wp*(x-xp(j))/((x1(i)-xp(j)));
        end
        W(x)=W+y1(i)*Wp;
    end
    W(x)=expand(W);
end

%% Literatura
% [1] Pańczyk B. i in.; Metody Numeryczne w przykładach, Politechnika Lubelska, Lublin 2012
% [2] Fortuna Z., Macukow B., Wąsowski J.; Metody Numeryczne; Wydawnictwa Naukowo-Techniczne, Warszawa
% 2001