clear all
close all
%
% compare quadrature formulas for the integral
%
% int_0^1 f(x)\phi_i(x)dx, i=2,m-1
%
mrange = 20:20:100;
order = [3/2,5/2,3,5];
p = [1,2,3,5];
for k = 1:4
  f = @(x) sqrt(abs(x-1/2)).^(2*p(k)-1); % in C^(p(k)-1)
  counter = 0;
  for m = mrange+1
    counter = counter+1;
    x = linspace(0,1,m)';
%    x = (-cos(linspace(0,pi,m)')+1)/2;
    h = diff(x);
% mass matrix    
    M = spdiags([[h/6;0],[h(1)/3;h(1:m-2)/3+h(2:m-1)/3;h(m-1)/3],[0;h/6]],[-1,0,1],m,m);
% trapezoidal
    qtrap = (f(x(1:m-2)) * 0 + f(x(2:m-1)) * 1) .* h(1:m-2) / 2 + ...
	(f(x(2:m-1)) * 1 + f(x(3:m)) * 0) .* h(2:m-1) / 2;
% barycentric
    qbar = (f(x(1:m-2))+f(x(2:m-1)))/2.*h(1:m-2)/2+...
	(f(x(2:m-1))+f(x(3:m)))/2.*h(2:m-1)/2;
% interpolation
    qinterp = M(2:m-1,:)*f(x);
% Gauss-Legendre
    xi1 = x(1:m-2)+(-sqrt(1/3)+1)/2*(x(2:m-1)-x(1:m-2));
    xi2 = x(1:m-2)+(sqrt(1/3)+1)/2*(x(2:m-1)-x(1:m-2));
    qGL = (f(xi1).*(xi1-x(1:m-2))/h+f(xi2).*(xi2-x(1:m-2))/h)*h/2;
    xi1 = x(2:m-1)+(-sqrt(1/3)+1)/2*(x(3:m)-x(2:m-1));
    xi2 = x(2:m-1)+(sqrt(1/3)+1)/2*(x(3:m)-x(2:m-1));	
    qGL = qGL+(f(xi1).*(x(3:m)-xi1)/h+f(xi2).*(x(3:m)-xi2)/h)*h/2;
    qexact = zeros(m-2,1);
    for i = 2:m-1
      integrand = @(xx) f(xx).*(xx-x(i-1))/h(i-1);
      qexact(i-1) = quadl(integrand,x(i-1),x(i),1e-14);
      integrand = @(xx) f(xx).*(x(i+1)-xx)/h(i);
      qexact(i-1) = qexact(i-1)+quadl(integrand,x(i),x(i+1),1e-14);
    end
      errtrap(counter) = norm(qexact-qtrap,inf);
      errbar(counter) = norm(qexact-qbar,inf);
      errinterp(counter) = norm(qexact-qinterp,inf);
      errGL(counter) = norm(qexact-qGL,inf);
  end
  figure
  loglog(mrange,mrange.^(-order(k)),'k',mrange,errbar,'*g',...
	 mrange,errtrap,'xr',mrange,errinterp,'ob',mrange,errGL,'+m')
  axis([10,1e3,1e-11,1e-1])
  legend(sprintf('order %.1f',order(k)),'barycentric',...
	 'trapezoidal','interpolation','GL2')
  xlabel('m')
  ylabel('error')
  print('-depsc',sprintf('quadrature1d%d.eps',k))
end
