clear all
close all
%
% compare quadrature formulas for the integral
%
% int_0^1 g(x)\phi_i(x)dx, i=2,m-1
%
p = 3; % smooth case g\in C^2
g = @(x) sqrt(abs(x-1/2)).^(2*p-1);
mrange = 50:100:450;
counter = 0;
for m = mrange+1
  counter = counter+1;
  x = linspace(0,1,m)';
  h = x(2)-x(1);
% trapezoidal
  qtrap = (g(x(1:m-2)) * 0 + g(x(2:m-1)) * 1) * h / 2 + ...
      (g(x(2:m-1)) * 1 + g(x(3:m)) * 0) * h / 2;
% barycentric
  qbar = (g(x(1:m-2))+g(x(2:m-1)))/2*h/2+...
      (g(x(2:m-1))+g(x(3:m)))/2*h/2;
  qexact = zeros(m-2,1);
  for i = 2:m-1
    integranda = @(xx) g(xx).*(xx-x(i-1))/h;
    qexact(i-1) = quadl(integranda,x(i-1),x(i),1e-12);
    integranda = @(xx) g(xx).*(x(i+1)-xx)/h;
    qexact(i-1) = qexact(i-1)+quadl(integranda,x(i),x(i+1),1e-12);
  end
  errtrap(counter) = norm(qexact-qtrap,inf);
  errbar(counter) = norm(qexact-qbar,inf);
end
loglog(mrange,mrange.^(-3),mrange,errbar,'*g',mrange,errtrap,'xr')
axis([10,1e3,1e-9,1e-2])
legend('order 3','barycentric','trapezoidal')
xlabel('m')
ylabel('error')
print('-depsc','qbarsmooth.eps')
p = 1; % less smooth case g\in C^0
g = @(x) sqrt(abs(x-1/2)).^(2*p-1);
counter = 0;
for m = mrange+1
  counter = counter+1;
  x = linspace(0,1,m)';
  h = x(2)-x(1);
% trapezi
  qtrap = 2 * (g(x(2:m-1)) * h / 2);
% bar
  qbar = (g(x(1:m-2))+g(x(2:m-1)))/2*h/2+...
      (g(x(2:m-1))+g(x(3:m)))/2*h/2;
  qexact = zeros(m-2,1);
  for i = 2:m-1
    integranda = @(xx) g(xx).*(xx-x(i-1))/h;
    qexact(i-1) = quadl(integranda,x(i-1),x(i));
    integranda = @(xx) g(xx).*(x(i+1)-xx)/h;
    qexact(i-1) = qexact(i-1)+quadl(integranda,x(i),x(i+1));
  end
  errtrap(counter) = norm(qexact-qtrap,inf);
  errbar(counter) = norm(qexact-qbar,inf);
end
figure
loglog(mrange,mrange.^(-3/2),mrange,errtrap,'xr',mrange,errbar,'*g')
axis([10,1e3,1e-9,1e-2])
legend('order 3/2','trapezoidal','barycentric')
xlabel('m')
ylabel('error')
print('-depsc','qbarnonsmooth.eps')