% Section 4.3 Cubic splines

clear

% sin(x) example with 4 nodes
x=[0 pi/2 pi 3*pi/2]';

A=zeros(4*3,4*3); y=zeros(4*3,1);
A(1,1:4)=x(1).^(3:-1:0);    y(1)=sin(x(1));
A(2,1:4)=x(2).^(3:-1:0);    y(2)=sin(x(2));
A(3,5:8)=x(2).^(3:-1:0);    y(3)=sin(x(2));
A(4,5:8)=x(3).^(3:-1:0);    y(4)=sin(x(3));
A(5,9:12)=x(3).^(3:-1:0);   y(5)=sin(x(3));
A(6,9:12)=x(4).^(3:-1:0);   y(6)=sin(x(4));
A(7,1:3)=[3*x(2)^2 2*x(2) 1];A(7,5:7)=-[3*x(2)^2 2*x(2) 1];  % y(7)=0;
A(8,5:7)=[3*x(3)^2 2*x(3) 1];A(8,9:11)=-[3*x(3)^2 2*x(3) 1]; % y(8)=0;
A(9,1:2)=[6*x(2) 2];A(9,5:6)=-[6*x(2) 2];   % y(9)=0;
A(10,5:6)=[6*x(3) 2];A(10,9:10)=-[6*x(3) 2]; % y(10)=0;

% natural BC
A(11,1:2)=[6*x(1) 2];    y(11)=0;
A(12,9:10)=[6*x(end) 2]; y(12)=0;
a=A\y;

t0=(0:.01:pi/2)';
p0=polyval(a(1:4),t0);
t1=(pi/2:.01:pi)';
p1=polyval(a(5:8),t1);
t2=(pi:.01:3*pi/2)';
p2=polyval(a(9:12),t2);
t=[t0;t1;t2];
figure;plot(t0,p0,'b',t1,p1,'b',t2,p2,'b',t,sin(t),'k',x,sin(x),'ko')

% clamped BC
A(11:12,:)=0; % erase previous BC
A(11,1:3)=[3*x(1)^2 2*x(1) 1];   y(11)=cos(x(1));
A(12,9:11)=[3*x(end)^2 2*x(end) 1];  y(12)=cos(x(end));
a=A\y;

p0=polyval(a(1:4),t0);
p1=polyval(a(5:8),t1);
p2=polyval(a(9:12),t2);
hold on;plot(t0,p0,'r',t1,p1,'r',t2,p2,'r')

% not-a-knot BC
A(11:12,:)=0; % erase previous BC
A(11,1)=1;A(11,5)=-1;    y(11)=0;
A(12,5)=1;A(12,9)=-1;    y(12)=0;
a=A\y;

p0=polyval(a(1:4),t0);
p1=polyval(a(5:8),t1);
p2=polyval(a(9:12),t2);
hold on;plot(t0,p0,'m-.',t1,p1,'m-.',t2,p2,'m-.')

% built-in spline
s=spline(x,sin(x),t);
figure;plot(t,s,'m',t,sin(t),'k',x,sin(x),'ko');