
%% opracowala Ewa Pabisek

function abc()
clear
format compact

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

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

%  definicja funkcji bazowych Hermite'a

 H1= inline('1 - 3*ksi.^2 + 2*ksi.^3');
 H2= inline('le*(ksi - 2*ksi.^2 + ksi.^3)', 'ksi','le');
 H3= inline('3*ksi.^2 - 2*ksi.^3','ksi');
 H4= inline('le*(ksi.^3 - ksi.^2)','ksi','le');


 a= 0; b= pi;   % rozwazany przedzial
 Lel= 4;        % liczba elementow 
 Lwez= Lel+1;  % liczba wezlow
 M= 9;          % przykladowa liczba punktow w elemencie


%% 2 %% GENEROWANIE WSPOLRZEDNYCH WEZLOW I WARTOSCI FUNKCJI    
 x= linspace(a,b,Lwez);
 y= F(x);  yp= Fp(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 Hermite''a')

%% 4 %%  ROZWIAZANIE ZADANIA
%       rysunek funkcji interpolujacych w podprzedzialach 
%       rysunek funkcji bledu
%       rysunek funkcji bazowych Hermite'a
 errnorm=0;
 Fnorm=0;
 for iel=1:Lel
    x1= x(iel); x2= x(iel+1);
    y1= y(iel); yp1= yp(iel); y2= y(iel+1); yp2= yp(iel+1);
    Le= x2-x1;

    xel= linspace(x1, x2, M);   
    ksi=(xel-x1)/Le;
    h1= H1(ksi) ;    h2= H2(ksi, Le) ;          % funkcje ksztaltu
    h3= H3(ksi);     h4= H4(ksi, Le) ; 
    yinter= y1*h1 +  yp1*h2 +  y2*h3 + yp2*h4;  % interpolant
    plot(xel, yinter, 'k-', 'LineWidth', 2)

    error= F(xel) - yinter;
    plot(xel, error, 'r', 'linewidth', 2);

    plot(xel, h1,'m', xel, h2,'g', xel, h3,'m--', xel, h4,'g--')
    
%  calkowanie metoda trapezow    
    hh= (x2-x1)/(M-1);
    error= error.^2;
    errnorm= errnorm + (sum(error)- error(1)/2 - error(M)/2) * hh;
    Fnorm= Fnorm + (sum(F(xel).^2) - F(x1)^2 - F(x2)^2) * hh;

    pause(1)
 end
 grid on

 errnorm
 Fnorm
 Rel_Err = sqrt(errnorm)/sqrt(Fnorm)*100%
 legend('Elem, Wez','F(x)','W_H(x_e)','Error','Location','southwest')
end



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

