
%% opracowala Ewa Pabisek

function abc()
 clear
 format compact

% 1 %%% DANE
     % zdefiniowac interpolowana funkcje F {inline('przepis')}

 F= inline('sin(3*x).*cos(x)','x')

 a= 0; b= pi;   % przedzial            
 Lel= 4;        % liczba elementów  
 Lwez= Lel+1;   % liczba wezlow   
 M= 9;          % liczba punktów w ES                


%% 2 %%% GENEROWANIE WSPLRZEDNYCH WEZLOW I WARTOSCI FUNKCJI  
 x= linspace(a,b,Lwez);
 y= F(x);  

 figure(1)
 clf
 plot(x,zeros(Lwez,1), 'k.-', 'LineWidth', 2, 'markersize', 30);  
 hold on, grid on; % axis equal,
 axis([a-.1 b+.1  -1 1.1]);

%% 3 %%% RYSUNEK FUNKCJI INTERPOLOWANEJ
 s= a : 0.01 : b;
 t= F(s);
 plot(s, t, '-c', 'LineWidth',3);  %,hold on, axis equal, grid on
 title(' Interpolacja Lagrange''a');

%% 4 %%% ROZWIAZANIE ZADANIA
%       rysunek funkcji interpolujacych w podprzedzialach 
%       rysunek funkcji bledu
%       rysunek funkcji bazowych Lagrange'a
 errnorm=0;
 Fnorm=0;
 for iel=1:Lel
    xL= x(iel); xP= x(iel+1);
    yL= y(iel); yP= y(iel+1);
    
    xel= linspace(xP, xL, M);   
    psi1= (xel-xP)/(xL-xP);     % funkcje ksztaltu
    psi2= (xel-xL)/(xP-xL);
    yinter= yL*psi1 + yP*psi2;  % interpolant
    plot(xel, yinter, 'k-', 'LineWidth', 2)

    error= F(xel) - yinter;
    plot(xel, error, 'r', 'LineWidth', 2);
    plot(xel, psi1,  'm', xel, psi2,  'g')
    
% kwadrat normy L2 bledu  i danej funkcji
% calkowanie metoda trapezów 
    hh= (xP-xL)/(M-1);
    error= error.^2;
    errnorm= errnorm + (sum(error)- error(1)/2 - error(M)/2) * hh;
    Fnorm= Fnorm + (sum(F(xel).^2) - F(xL)^2 - F(xP)^2) * hh;

    pause(1)
 end
 legend('ElemWez','F(x)','W_L(x_e)','Error','Location','southwest')

 Rel_Err = sqrt(errnorm)/sqrt(Fnorm)*100%

end


%e(x)_2= sqrt{ calka_x1^x2  [F(x) - L(x)]^2 dx }
%calka_tr= H [1/2 f_0 + sum_1^n-1 f_j + 1/2 f_N]

