function fig5_19

% error when fitting with cubic splines and cheb points
% integration rule is composite trap (see first edition, exercise 6.25)

%  determine exact values
a=-1; b=1;
Ne=8000;
xp=linspace(a,b,Ne);
h=xp(2)-xp(1);
for i=1:Ne
    R(i)=f(xp(i));
end

N=7;

tic
%  spline
nx=5;
for in=1:N
    nx=2*nx;
    num(in)=nx;
    xd=linspace(a,b,nx);
    for iy=1:nx
        yd(iy)=f(xd(iy));
    end
    xp=linspace(a,b,Ne);
    % natural spline
    ys = splineA(xd,yd,xp);
    eS(in)=h*norm(ys-R,1);
    % clamped spline
    yf = splineA(xd,yd,xp,[1 1],[50/26^2 -50/26^2]);
    ef(in)=h*norm(yf-R,1);
end
toc

tic
%  cheb
nx=5;
for in=1:N
    nx=2*nx;
    num(in)=nx;
    xd=linspace(a,b,nx);
    for ix=1:nx
        cc=(2*ix-1)*pi/(2*nx);
        xxd(ix)=0.5*(a+b+(b-a)*cos(cc));
        yyd(ix)=f(xxd(ix));
    end
    data=[xxd' yyd'];
    [yc]=barylag(data,xp');
    eC(in)=h*(norm(yc'-R,1)-0.5*abs(yc(1)-R(1))-0.5*abs(yc(Ne)-R(Ne)));
end
toc

clf
% get(gcf)
set(gcf,'Position', [25 1115 658 230])
subaxis(1,1,1,1,'MT',0.01,'MB',0.19,'MR',-0.01,'ML',0.07,'P',0.04)
co = [0 0 1;
    0 0.5 0;
    1 0 0;
    0 0.75 0.75;
    0.75 0 0.75;
    0.75 0.75 0;
    0.25 0.25 0.25];
set(groot,'defaultAxesColorOrder',co)

loglog(num,eS,'--or','LineWidth',1.6,'MarkerSize',8)
hold on
loglog(num,eC,'--*b','LineWidth',1.6,'MarkerSize',8)
%loglog(num,ef,'--dk','LineWidth',1.6,'MarkerSize',8)
box on
grid on
axis([10 1e3 1e-16 5e-1])
set(gca,'YTick',[1e-16 1e-11 1e-6 1e-1 1])
xlabel('n-axis')
ylabel('Error')

legend({' Spline',' Cheb'},'Location','SouthWest','FontSize',16,'FontWeight','bold')
%legend({' Spline',' Cheb',' Fixed'},'Location','SouthWest','FontSize',16,'FontWeight','bold')

set(gca,'FontSize',16,'FontWeight','bold')

%exportgraphics(gcf,'/Users/mark/Desktop/splineE.eps')

function g=f(x)
g=1/(1+25*x^2);

function p=ell(i,x,xd)
[n1 n2]=size(xd);
p=1;
for j=1:n2
    if j ~= i
        p=p*(x-xd(j))/(xd(i)-xd(j));
    end
end


function [p]=barylag(data,x)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% barylag.m
%
% Interpolates the given data using the Barycentric
% Lagrange Interpolation formula. Vectorized to remove all loops
%
% data - a two column vector where column one contains the
%        nodes and column two contains the function value
%        at the nodes
% p - interpolated data. Column one is just the
%     fine mesh x, and column two is interpolated data
%
% Reference:
%
% (1) Jean-Paul Berrut & Lloyd N. Trefethen, "Barycentric Lagrange
%     Interpolation"
%     http://web.comlab.ox.ac.uk/oucl/work/nick.trefethen/berrut.ps.gz
% (2) Walter Gaustschi, "Numerical Analysis, An Introduction" (1997) pp. 94-95
%
%
% Written by: Greg von Winckel       03/07/04
% Contact:    gregvw@chtm.unm.edu
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

M=length(data);     N=length(x);

% Compute the barycentric weights
X=repmat(data(:,1),1,M);

% matrix of weights
W=repmat(1./prod(X-X.'+eye(M),1),N,1);

% Get distances between nodes and interpolation points
xdist=repmat(x,1,M)-repmat(data(:,1).',N,1);

% Find all of the elements where the interpolation point is on a node
[fixi,fixj]=find(xdist==0);

% Use NaNs as a place-holder
xdist(fixi,fixj)=NaN;
H=W./xdist;

% Compute the interpolated polynomial
p=(H*data(:,2))./sum(H,2);

% Replace NaNs with the given exact values.
p(fixi)=data(fixj,2);




