% Math 5610 - HW 5 - Problem 1 - SAMPLE CODE
%
% Spline interpolation - B&F 3.5.3 d and 3.5.5 d
%
% WHAT YOU NEED TO DO:
%  * change x,y, f and xeval to data in 3.5.4 d and 3.5.6 d
%  * figure out interval index j where evaluation point xeval lies
%  * how to evaluate derivative of spline at eval point
%    xeval (which happens to be 0.25 for both)

% nodes and values
n = 3;
x = [ 0.1, 0.2, 0.3, 0.4];
y = [ -0.62049958, -0.28398668, 0.00660095, 0.24842440];

% evaluation node
xeval=0.25;

% h(i) contains x(i+1)-x(i)
h = diff(x);

% underlying function and derivative
f = @(x) x.*cos(x) -2*x.^2 + 3*x -1;
fp = @(x) cos(x) -x.*sin(x) -4*x + 3;

% coefficients a
a = y;
% form linear system from theorem 3.11
A = [ 1, 0, 0, 0
      h(1), 2*(h(1)+h(2)),  h(2), 0
      0, h(2), 2*(h(2)+h(3)),  h(3)
      0, 0, 0, 1];
B = [0
     (3/h(2))*(a(3)-a(2)) - (3/h(1))*(a(2)-a(1))
     (3/h(3))*(a(4)-a(3)) - (3/h(2))*(a(3)-a(2))
     0];

% compute coefficients
c = A\B;

% get coefficients b and d
b = zeros(n,1); d = zeros(n,1);
for j=1:n,
 b(j) = (a(j+1)-a(j))/h(j) - h(j)*(c(j+1)+2*c(j))/3;
 d(j) = (c(j+1)-c(j))/(3*h(j));
end;

%%%%%%%%%% Table with spline coefficients %%%%%%%%%%%%%%
fprintf('%2s %13s %13s %13s %13s\n','i','a','b','c','d');
for j=1:n,
 fprintf('%2d %13.8g %13.8g %13.8g %13.8g\n',j-1,a(j),b(j),c(j),d(j));
end;

%%%%%%%%%% evaluate spline at xeval and compare to function %%%%%%%%%%%%%%
j = 2; % interval index where xeval lies
yy = a(j) + b(j)*(xeval-x(j)) + c(j)*(xeval-x(j)).^2 + d(j)*(xeval-x(j)).^3;
fprintf('%5s %13s %13s %13s\n','x','S(x)','f(x)','|S(x)-f(x)|');
fprintf('%5.2f %13.8g %13.8g %13.4e\n',xeval,yy,f(xeval),abs(yy-f(xeval)));

%%%%%%%%%% evaluate spline derivative at xeval and compare to function %%%%%%%%%%%%%%
yy =  TO BE COMPLETED
fprintf('%5s %13s %13s %13s\n','x','S''(x)','f''(x)','|S''(x)-f''(x)|');
fprintf('%5.2f %13.8g %13.8g %13.4e\n',xeval,yy,fp(xeval),abs(yy-fp(xeval)));

%%%%%%%%%% evaluate and plot spline %%%%%%%%%%%%%%
figure(1); clf; hold on;
% plot data nodes
plot(x,y,'k*');
nf = 10; % number of points used to plot in a subinterval
for j=1:n,
 xx = linspace(x(j),x(j+1),nf);
 yy = a(j) + b(j)*(xx-x(j)) + c(j)*(xx-x(j)).^2 + d(j)*(xx-x(j)).^3;
 plot(xx,yy,'b-');
end;
hold off;
xlabel('x'); ylabel('y'); title('S(x)');

% Sample  output
%   >> prob1
%    i             a             b             c             d
%    0   -0.62049958     3.4550869             0    -8.9957933
%    1   -0.28398668     3.1852131     -2.698738   -0.94630333
%    2    0.00660095     2.6170764     -2.982629     9.9420967
%       x          S(x)          f(x)   |S(x)-f(x)|
%    0.25   -0.13159116   -0.13277189    1.1807e-03
%       x         S'(x)         f'(x) |S'(x)-f'(x)|
%    0.25     2.9082421     2.9070614    1.1806e-03
%   >> 

