clear all
close all
mrange = 64:100:664;
a = 0;
b = 2*pi;
xx = linspace(a,b,6000);
f = @(x) cos(x)+(x-2*pi).^2/2;
f1 = @(x) -sin(x)+(x-2*pi);
f2 = @(x) -cos(x)+1;
f3 = @(x) sin(x);
%f = @(x) abs(x).^(7/2);
%f1 = @(x) 7/2*abs(x).^(3/2).*x;
%f2 = @(x) 35/4*abs(x).^(3/2);
%f3 = @(x) 105/8*abs(x).^(1/2).*sign(x);
% f(x)
counter = 0;
for m = mrange;
  counter = counter+1;
  x = linspace(a,b,m);
  h(counter) = (b-a)/(m-1);
  y = f(x);
  pp = csape(x,y,'variational');
  yy = ppval(pp,xx);
  errnat(counter) = norm(yy-f(xx),inf);
  pp = csape(x,y,'not-a-knot');
  yy = ppval(pp,xx);
  errnak(counter) = norm(yy-f(xx),inf);
  pp = csape(x,y,'periodic');
  yy = ppval(pp,xx);
  errper(counter) = norm(yy-f(xx),inf);
  pp = csape(x,[f1(a),y,f1(b)],'complete');
  yy = ppval(pp,xx);
  errvin(counter) = norm(yy-f(xx),inf);
end
loglog(h,errnat,'*',h,errnak,'o',...
       h,errper,'x',h,errvin,'+')
grid
axis([1e-3,1,1e-12,1e4])
legend('natural','not-a-knot',...
'periodic','vincolated','location','northwest')
xlabel('h')
ylabel('error in infinity norm')
title('f(x)')
% f'(x)
counter = 0;
for m = mrange;
  counter = counter+1;
  x = linspace(a,b,m);
  h(counter) = (b-a)/(m-1);
  y = f(x);
  pp = csape(x,y,'variational');
  [breaks,coefs] = unmkpp(pp);
  coefs = [3*coefs(:,1),2*coefs(:,2),coefs(:,3)];
  pp = mkpp(breaks,coefs);
  yy = ppval(pp,xx);
  errnat(counter) = norm(yy-f1(xx),inf);
  pp = csape(x,y,'not-a-knot');
  [breaks,coefs] = unmkpp(pp);
  coefs = [3*coefs(:,1),2*coefs(:,2),coefs(:,3)];
  pp = mkpp(breaks,coefs);
  yy = ppval(pp,xx);
  errnak(counter) = norm(yy-f1(xx),inf);
  pp = csape(x,y,'periodic');
  [breaks,coefs] = unmkpp(pp);
  coefs = [3*coefs(:,1),2*coefs(:,2),coefs(:,3)];
  pp = mkpp(breaks,coefs);
  yy = ppval(pp,xx);
  errper(counter) = norm(yy-f1(xx),inf);
  pp = csape(x,[f1(a),y,f1(b)],'complete');
  [breaks,coefs] = unmkpp(pp);
  coefs = [3*coefs(:,1),2*coefs(:,2),coefs(:,3)];
  pp = mkpp(breaks,coefs);
  yy = ppval(pp,xx);
  errvin(counter) = norm(yy-f1(xx),inf);
end
figure
loglog(h,errnat,'*',h,errnak,'o',...
       h,errper,'x',h,errvin,'+')
grid
axis([1e-3,1,1e-12,1e4])
legend('natural','not-a-knot',...
'periodic','vincolated','location','northwest')
xlabel('h')
ylabel('error in infinity norm')
title('f''(x)')
% f''(x)
counter = 0;
for m = mrange;
  counter = counter+1;
  x = linspace(a,b,m);
  h(counter) = (b-a)/(m-1);
  y = f(x);
  pp = csape(x,y,'variational');
  [breaks,coefs] = unmkpp(pp);
  coefs = [6*coefs(:,1),2*coefs(:,2)];
  pp = mkpp(breaks,coefs);
  yy = ppval(pp,xx);
  errnat(counter) = norm(yy-f2(xx),inf);
  pp = csape(x,y,'not-a-knot');
  [breaks,coefs] = unmkpp(pp);
  coefs = [6*coefs(:,1),2*coefs(:,2)];
  pp = mkpp(breaks,coefs);
  yy = ppval(pp,xx);
  errnak(counter) = norm(yy-f2(xx),inf);
  pp = csape(x,y,'periodic');
  [breaks,coefs] = unmkpp(pp);
  coefs = [6*coefs(:,1),2*coefs(:,2)];
  pp = mkpp(breaks,coefs);
  yy = ppval(pp,xx);
  errper(counter) = norm(yy-f2(xx),inf);
  pp = csape(x,[f1(a),y,f1(b)],'complete');
  [breaks,coefs] = unmkpp(pp);
  coefs = [6*coefs(:,1),2*coefs(:,2)];
  pp = mkpp(breaks,coefs);
  yy = ppval(pp,xx);
  errvin(counter) = norm(yy-f2(xx),inf);
end
figure
loglog(h,errnat,'*',h,errnak,'o',...
       h,errper,'x',h,errvin,'+')
grid
axis([1e-3,1,1e-12,1e4])
legend('natural','not-a-knot',...
'periodic','vincolated','location','northwest')
xlabel('h')
ylabel('error in infinity norm')
title('f''''(x)')
% f'''(x)
counter = 0;
for m = mrange;
  counter = counter+1;
  x = linspace(a,b,m);
  h(counter) = (b-a)/(m-1);
  y = f(x);
  pp = csape(x,y,'variational');
  [breaks,coefs] = unmkpp(pp);
  coefs = [6*coefs(:,1)];
  pp = mkpp(breaks,coefs);
  yy = ppval(pp,xx);
  errnat(counter) = norm(yy-f3(xx),inf);
  pp = csape(x,y,'not-a-knot');
  [breaks,coefs] = unmkpp(pp);
  coefs = [6*coefs(:,1)];
  pp = mkpp(breaks,coefs);
  yy = ppval(pp,xx);
  errnak(counter) = norm(yy-f3(xx),inf);
  pp = csape(x,y,'periodic');
  [breaks,coefs] = unmkpp(pp);
  coefs = [6*coefs(:,1)];
  pp = mkpp(breaks,coefs);
  yy = ppval(pp,xx);
  errper(counter) = norm(yy-f3(xx),inf);
  pp = csape(x,[f1(a),y,f1(b)],'complete');
  [breaks,coefs] = unmkpp(pp);
  coefs = [6*coefs(:,1)];
  pp = mkpp(breaks,coefs);
  yy = ppval(pp,xx);
  errvin(counter) = norm(yy-f3(xx),inf);
end
figure
loglog(h,errnat,'*',h,errnak,'o',...
       h,errper,'x',h,errvin,'+')
grid
axis([1e-3,1,1e-12,1e4])
legend('natural','not-a-knot',...
'periodic','vincolated','location','northwest')
xlabel('h')
ylabel('error in infinity norm')
title('f''''''(x)')
