% Beispiel Polynominterpolation <> kubischer Spline
% M.Ries, Mai 2019
% Systemmatrix
A = [1 0 0 0  0  0  0  0 ;  % s_0(0)  = -1 , Startwert s_0
     0 0 0 0  1  1  1  1 ;  % s_1(1)  =  0 , Startwert s_1
     0 0 0 0  1  2  4  8 ;  % s_1(2)  =  0 , Endwert s_1
  %      zusammen m+1
     0 0 2 0  0  0  0  0 ;  % S"_0(0) = 0  , natürliche Bedingung
     0 0 0 0  0  0  2 12 ;  % S"_1(2) = 0  ,   =    =   =
  %      also 2
     1 1 1 1 -1 -1 -1 -1 ;  % s_0(1)  = s_1(1) , Funktionswert
  %      zusammen m-1
     0 1 2 3  0 -1 -2 -3 ;  % s'_0(1) = s'_1(1) , Tangente
  %      nochmal m-1
     0 0 2 6  0  0 -2 -6 ]  % s"_0(1) = s"_1(1) , Krümmung
  %      und m-1, also insgesamt (m+1) + 2 + 3*(m-1) = 4m 
  %
% rechte Seite für die beiden Fälle
b_1 = [-1 0 1    0 0   0 0 0 ]'
b_2 = [-1 0 0    0 0   0 0 0 ]'

x_1 = mldivide( A, b_1)   % Lösung 
x_2 = mldivide( A, b_2)

t_0 = linspace( 0, 1) ;
t_1 = linspace( 1, 2) ;
% Die beiden Teil-Splinefunktionen

s_0 = flipud( x_2(1:4) )  % up down
ss_0 = polyval(s_0, t_0) ;
s_1 = flipud( x_2(5:8) )   % up down
ss_1 = polyval(s_1, t_1) ;
%%s_0 = -1.0     +  5.0/4.0*t_0                  - 1.0/4.0*t_0.^3 ;
%%s_1 = -3.0/2.0 + 11.0/4.0*t_1 - 3.0/2.0*t_1.^2 + 1.0/4.0*t_1.^3 ;

% Vergleich mit dem Lagrange-Polynom
p = 1.0/2.0 * fliplr([-2, 3, -1])  % left right
tt = [t_0, t_1] ;  % Das gesamte Intervall
pp = polyval(p, tt ) ;
%%p  = 1.0/2.0 * ( -tt.^2 + 3*tt - 2) ;

t_werte = [ 0  1  2 ]' ;
y_werte = [-1  0  0 ]' ;

%clf

figure('Name','Meßwerte','NumberTitle','off') 
hold on
axis([-0.5 2.5 -1.5 0.5])
plt_punkte = plot( t_werte, y_werte,'o') ; 
plt_lagrange = plot( tt, pp, 'r') ;   % Polynom

plt_s_0 = plot( t_0, ss_0, 'b') ; % s_0
plt_s_1 = plot( t_1, ss_1, 'b') ; % s_1

legend([plt_punkte plt_lagrange plt_s_0],{'Meßwerte','Lagrange','Spline'},'Location','southeast')