% Script 1 for Week 8

% Given independent data

t=[0 0.1 0.3 0.5 0.7 0.8 0.9 1.1 1.3 1.4 1.5 1.6 1.8 1.9];

% integrate the dependent data

f=@(x) sin(x)+cos(x)
z=f(t);

% over the interval.

% Of course we can get an exact answer 

format SHORTE

exact=integral(@(x) sin(x)+cos(x),0,1.9)

% exact =    2.2696

% which is the same as

-cos(1.9)+sin(1.9)+1

tr=trapz(t,z)

% tr = 2.2640

% with error

error(1)=(tr-exact)/exact

% ans =   -0.0025

% which is pretty good.

% Rather than using trapz we could create a spline interpolation function on
% this data and integrate THAT.

fspl=@(w) spline(t,z,w)

x=linspace(0,1.9);
yex=f(x);
yspl1=fspl(x);
plot(t,z,'*',x,yex,x,yspl1)


intfspl=integral(fspl,0,1.9)

% intfspl=    2.2696

error(2)=(intfspl-exact)/exact

% ans =  -1.8949e-06

% There are various options in the spline function-creator.

fclamp=@(w)spline(t,[-5 z 7],w);
yclmp=fclamp(x);

% This forces the spline to have derivative -5 at the left edge and 7 at 
% the right. Obviously not optimal for our situation. This is called 
% "clamping" the spline.  For example, the choice [0 z 0] enforces 
% zero derivatives at the endpoints.

plot(t,z,'*',x,yex,x,yclmp)

intfclmp=integral(fclamp,0,1.9)

% intfclmp =    2.2714

error(3)=(intfclmp-exact)/exact

% ans =   8.1077e-04

% Interpolation with the pchip (piecewise cubic Hermite) option may 
% be best for some purposes: it has less tendency to "overshoot" 
% but discontinuity of the second derivative might be visible.

fpchip=@(w)pchip(t,z,w);
ychip=fpchip(x);

plot(t,z,'*',x,yex,x,ychip)

intpchp=integral(fpchip,0,1.9)

% intpchp =    2.2695

error(4)=(intpchp-exact)/exact

% ans =  -2.1291e-05

% The natural spline interpolation uses, as the final two conditions to 
% define the coefficients, setting the 2nd derivative at the first and 
% last knots equal to zero. That means you could extend the spline 
% beyond the knots by a straight line and it would "look right" as far 
% as concavity "violations" are concerned.
% Use the posted natspline.m file for this. It seems that natspline.m 
% will not allow you to extrapolate: it insists that evaluations be 
% among the knots. However it will provide the derivative at the 
% endpoints, so you could extend linearly with this slope.

[a,b,c]=natspline(t,z,1.9)

% a = 0.6230   b = -1.2508  c = 0

% So the slope of the tangent at (1.9, 0.62301)  is  -1.2508 and the 2nd 
% derivative is 0 so we could extrapolate by extending linearly beyond 
% the data interval to the right using this slope, with continuous 2nd 
% derivative there.

fnat=@(w) natspline(t,z,w)

ynat=fnat(x);

plot(t,z,'*',x,yex,x,ynat)

intnat=integral(fnat,0,1.9)

% intnat =2.2696

error(5)=(intnat-exact)/exact

% ans =3.4952e-07


A few of the cubic interpolation definitions:

After forcing the cubic pieces to go through the data points and 
continuity of 1st and 2nd derivatives at the knots there are still 
two free variables. These can be chosen in various ways.

Definition of four piecewise-cubic interpolation methods:

not-a-knot cubic spline
force continuity of the 3rd derivative at the second and the next-to-last
knot. This is the default interpolation of the spline MATLAB command

clamped cubic spline  
force 1st derivative to be a specified value at the endpoints

natural cubic spline   
force 2nd derivative to be 0 at first and last knots (so could extend 
linearly past these knots)

pchip    
piecewice cubic Hermite polynomial: a different method entirely
2nd derivatives may not be continuous which could, potentially, be visible. 
1st derivatives at the knots will not be the same as for splines. With this method the interpolated values do not tend to overshoot the data as 
much as can happen with splines.

spline and interp1 with the spline option are essentially the same. spline 
performs interpolation on rows of an input matrix, while 
interp1(x,y,xq,'spline') performs interpolation on columns of an input matrix.

There is a way to make MATLAB produce the actual cubic coefficients on each 
piece, though it is unclear why you would want to. Investigate the unmkpp 
function for this.