%% Numerik D-MAVT, RW (Nipp), SS06
%
% This is the M-file solution to the seventh exercise sheet.
% To run, save this file to your account and ensure that the
% "Current directory" in MATLAB is the same as that where the
% file is saved.
%
%
% Then to run the file type "s7" in the command window


function s7
disp(['Aufgabe 4a'])
format compact;
x=0:.5:1.5; f=[0 1 0 -1]; % data
myspline(x,f,.25),

disp(['Aufgabe 4b'])
format compact;
x=-5+2*([1:6]-1);	% data
f=1./(1+x.^2);
xx=-4:2:4;		% where to eval spline
y=myspline(x,f,xx),	% compute and evaluate spline
% output:
% y = 
%       -0.0285    0.3131    0.5685    0.3131   -0.0285
ym=spline(x,f,xx),	% now use matlab function spline
% output:
% ym = 
%       -0.0285    0.3131    0.5685    0.3131   -0.0285
% make a pic
xp=-5:0.1:5; fp=1./(1+xp.^2); 	% the function
yp=myspline(x,f,xp);		% the spline
plot(xp,fp,'-.',x,f,'o',xx,y,'.',xx,ym,'*',xp(1:length(xp)-1),yp,':');
legend('f(x)','data','Q_i(x)','MATLAB','Spline');




function yy = myspline(x, f, xx)

% spline (not-a-knot) through data (x_i,f_i), evaluated at xx
% Numerik D-MAVT, RW (Nipp), SS 06

% Serie 7, A4
% in:	x,f	data
%	xx	evaluate resulting spline at xx
% out:	yy	spline evaluated at xx
% auth: fab@math.ethz.ch

n = length(x);	% number of data points
h = diff(x)';	% stepsizes h_i
% build matrix of linear system
b = 1./h;
a = 2*(b(1:n-2)+b(2:n-1));
c = diff(f)'./(h.^2);
A = sparse(diag([b(1);a(1:n-2);b(n-1)]) + diag([a(1)/2;b(2:n-1)],1) + diag([b(1:n-2);a(n-2)/2],-1));
% same as
% A = spdiags([[b(1:n-2);.5*a(n-2);0],[b(1);a;b(n-1)],[0;.5*a(1);b(1:n-2)] ],[-1,0,1], n, n),

% build rhs of linear system	    
d = [2*c(1)+h(1)/(h(1)+h(2))*(c(1)+c(2)); ...
     3*(c(1:n-2) + c(2:n-1)); ...
     2*c(n-1)+h(n-1)/(h(n-1)+h(n-2))*(c(n-1)+c(n-2))];
fp = A\d;	% solve linear system
% evaluate at xx
yy=[];
int=[]; % find intervals
for i=1:length(xx) int=[int max(find(x(1:length(x)-1)<xx(i)))]; end
for k=1:length(int)
  i = int(k);
  t = (xx(k)-x(i))/h(i);		% local coordinate
  % evaluate polynomial Q_i(t) (Hermite)
  a = [f(i),f(i+1)];
  b = [h(i)*fp(i), a(2)-a(1), h(i)*fp(i+1)];
  c = diff(b);
  d = c(2)-c(1);
  yy = [yy a(1) + (b(2)+(c(1)+d*t)*(t-1))*t]; % Hermite
end
