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)
mm = 1001;
xx = linspace(0,1,mm)';
mrange = 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 m = mrange+1
    counter = counter+1;
    x = linspace(0,1,m)';
    h = diff(x);
    n = m-1;
    ell = zeros(n,2);
    a = zeros(n,2);
    b = zeros(n,2);
    delta = zeros(n,1);
    fell = zeros(n,1);
% topology
    for j = 1:n
      ell(j,1) = j;
      ell(j,2) = j+1;
      a(j,1) = x(ell(j,2));
      b(j,1) = -1;
      a(j,2) = -x(ell(j,1));
      b(j,2) = 1;
      delta(j) = x(ell(j,2))-x(ell(j,1));
      fell(j) = (ffun(x(ell(j,1)))+ffun(x(ell(j,2))))/2;
    end
    A = spalloc(m,m,3*m-2);
    f = zeros(m,1);
    fexact = zeros(m,1);
    for j = 1:n
      for k = 1:2
	A(ell(j,k),ell(j,k)) = A(ell(j,k),ell(j,k))+1/delta(j);
	integrand = @(xx) ffun(xx).*(a(j,k)+b(j,k)*xx)/delta(j);
	tildefjk = quadl(integrand,x(ell(j,1)),x(ell(j,2)),1e-14);
	fexact(ell(j,k)) = fexact(ell(j,k))+tildefjk;
	tildefjk = fell(j)*h(j)/2; % barycentric formula
	f(ell(j,k)) = f(ell(j,k))+tildefjk; 
	for i = k+1:2
	  A(ell(j,k),ell(j,i)) = A(ell(j,k),ell(j,i))-1/delta(j);
	  A(ell(j,i),ell(j,k)) = A(ell(j,k),ell(j,i)); % symmetry     
	end
      end
    end
    M = 2^52; % penalization method
% Dirichlet b.c.
    A(1,1) = M;
    fexact(1) = 0*M;
    f(1) = 0*M;
    A(m,m) = M;
    fexact(m) = 0*M;
    f(m) = 0*M;
% solution
    uhat = A\fexact;
    ubar = A\f;
% 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(mm-1);
    errL2bar(counter) = norm(uubar-uexact(xx))/sqrt(mm-1);
    norm(uhat-uexact(x),inf) % uhat is interpolant!
  end
  figure
  loglog(mrange,mrange.^(-2),mrange,errL2bar,'*g',mrange,errL2,'ok')
  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
