clear all
close all
% solve the problem
%
% -u''(x) = f(x) = sqrt(abs(x-1/2))^(2*p-1)
%
% + Homogeneous Dirichlet b.c.
%
% The exact solution is 
%
% u(x) = -4/((2*p+3)*(2*p+1))*sqrt(abs(x-1/2))^(2*p+3)-sqrt(1/2)^(2*p+3)
%
% f\in H^(p-1)\cap C^(p-1) % not so regular
% u\in H^(p+1)\cap C^(p+1)
nn = 1001;
xx = linspace(0,1,nn)';
nrange = 50:100:450;
p = [1,3]; 
for l = 1:2
  ffun = @(x) sqrt(abs(x-1/2)).^(2*p(l)-1);
  uexact = @(x) -4/((2*p(l)+3)*(2*p(l)+1))*...
      (sqrt(abs(x-1/2)).^(2*p(l)+3)-sqrt(1/2)^(2*p(l)+3));
  counter = 0;
  for n = nrange+1
    counter = counter+1;
    x = linspace(0,1,n)';
    h = diff(x);
    M = n-1;
    ell = zeros(M,2);
    alpha = zeros(M,2);
    beta = zeros(M,2);
    delta = zeros(M,1);
    fell = zeros(M,1);
% topology
    for m = 1:M
      ell(m,1) = m;
      ell(m,2) = m+1;
      alpha(m,1) = x(ell(m,2));
      beta(m,1) = -1;
      alpha(m,2) = -x(ell(m,1));
      beta(m,2) = 1;
      delta(m) = x(ell(m,2))-x(ell(m,1));
      fell(m) = (ffun(x(ell(m,1)))+ffun(x(ell(m,2))))/2;
    end
    a = spalloc(n,n,3*n-2);
    ftilde = zeros(n,1);
    f = zeros(n,1);
    for m = 1:M
      for k = 1:2
	a(ell(m,k),ell(m,k)) = a(ell(m,k),ell(m,k))+1/delta(m);
	integrand = @(xx) ffun(xx).*(alpha(m,k)+beta(m,k)*xx)/delta(m);
	f(ell(m,k)) = f(ell(m,k))+...
		      quadl(integrand,x(ell(m,1)),x(ell(m,2)),1e-14);
	tildeFmk = fell(m)*h(m)/2; % barycentric formula
	ftilde(ell(m,k)) = ftilde(ell(m,k))+tildeFmk; 
	for i = k+1:2
	  a(ell(m,k),ell(m,i)) = a(ell(m,k),ell(m,i))-1/delta(m);
	  a(ell(m,i),ell(m,k)) = a(ell(m,k),ell(m,i)); % symmetry     
	end
      end
    end
    tgv = 2^52; % penalization method
% Dirichlet b.c.
    a(1,1) = tgv;
    f(1) = 0*tgv;
    ftilde(1) = 0*tgv;
    a(n,n) = tgv;
    f(n) = 0*tgv;
    ftilde(n) = 0*tgv;
% solution
    uhat = a\f;
    ubar = a\ftilde;
% linear interpolation
    uuhat = interp1(x,uhat,xx,'linear');
    uubar = interp1(x,ubar,xx,'linear');
% approximation of the L^2 norm (by trapezoidal rule on xx)
    errL2(counter) = norm(uuhat-uexact(xx))/sqrt(nn-1);
    errL2bar(counter) = norm(uubar-uexact(xx))/sqrt(nn-1);
    norm(uhat-uexact(x),inf) % uhat is interpolant!
  end
  figure
  loglog(nrange,nrange.^(-2),nrange,errL2bar,'*g',nrange,errL2,'ok')
  title(sprintf('f in C^%d',p(l)-1))
  legend('order 2','error in L2 norm','error in L2 norm (exact rhs)')
  axis([10,1e3,1e-8,1e-3])
  xlabel('m')
  ylabel('error')
  print('-depsc',sprintf('assembly1d%d.eps',n))
end
