clear all
close all
% solve the problem
% -u''(x) = -(2*p+3)*(2*p+1)/4*sqrt(abs(x-1/2))^(2*p-1)
% + Homogeneous Dirichlet b.c.
%
% The exact solution is u(x) = sqrt(abs(x-1/2))^(2*p+3)-(1/2)^(2*p+3)
% u\in H^(p+1)
mm = 1001;
xx = linspace(0,1,mm)';
p = 3; % smooth case
uexact = @(x) sqrt(abs(x-1/2)).^(2*p+3)-sqrt(1/2)^(2*p+3);
mrange = 50:100:450;
counter = 0;
for m = mrange+1
 counter = counter+1;
 x = linspace(0,1,m)';
 h = diff(x);
 l = m-1;
 gfun = @(x) -(2*p+3)*(2*p+1)/4*sqrt(abs(x-1/2)).^(2*p-1);
 ell = zeros(l,2);
 gell = zeros(l,1);
 for j = 1:l
   ell(j,1) = j;
   ell(j,2) = j+1;
   gell(j) = (gfun(x(ell(j,1)))+gfun(x(ell(j,2))))/2;
 end
 A = spalloc(m,m,3*(m-2)+2*2);
 g = zeros(m,1);
 gexact = zeros(m,1);
 for j = 1:l
   for k = 1:2
     A(ell(j,k),ell(j,k)) = A(ell(j,k),ell(j,k))+1/h(j);
     integranda = @(xx) gfun(xx).*...
	 ((2-k)*(x(j+1)-xx)/h(j)+(k-1)*(xx-x(j))/h(j));
     gexact(ell(j,k)) = gexact(ell(j,k))+quadl(integranda,x(j),x(j+1),1e-12);
     g(ell(j,k)) = g(ell(j,k))+gell(j)*h(j)/2; % barycentric formula
     for i = k+1:2
       A(ell(j,k),ell(j,i)) = A(ell(j,k),ell(j,i))-1/h(j);
       A(ell(j,i),ell(j,k)) = A(ell(j,k),ell(j,i));      
     end
    end
  end
% Dirichlet b.c.
  A(1,1) = 1e16;
  gexact(1) = 0;
  g(1) = 0;
  A(m,m) = 1e16;
  gexact(m) = 0;
  g(m) = 0;
% solution
  u = A\gexact;
  ubar = A\g;
% linear interpolation
  uu = interp1(x,u,xx,'linear');
  uubar = interp1(x,ubar,xx,'linear');
% approximation of the L^2 norm (by trapezoidal rule on xx)
  errL2(counter) = norm(uu-uexact(xx))/sqrt(mm-1);
  errL2bar(counter) = norm(uubar-uexact(xx))/sqrt(mm-1);
%  errL2(counter) = norm(u-uexact(x))/sqrt(m-1); % interpolant!
%  errL2bar(counter) = norm(ubar-uexact(x))/sqrt(m-1);

end
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-7,1e-1])
xlabel('m')
ylabel('error')
print('-depsc','fem1dsmooth.eps')
p = 1; % less smooth case
uexact = @(x) sqrt(abs(x-1/2)).^(2*p+3)-sqrt(1/2)^(2*p+3);
counter = 0;
for m = mrange+1
 counter = counter+1;
 x = linspace(0,1,m)';
 h = diff(x);
 l = m-1;
 gfun = @(x) -(2*p+3)*(2*p+1)/4*sqrt(abs(x-1/2)).^(2*p-1);
 ell = zeros(l,2);
 gell = zeros(l,1);
 for j = 1:l
   ell(j,1) = j;
   ell(j,2) = j+1;
   gell(j) = (gfun(x(ell(j,1)))+gfun(x(ell(j,2))))/2;
 end
 A = spalloc(m,m,3*(m-2)+2*2);
 g = zeros(m,1);
 gexact = zeros(m,1);
 for j = 1:l
   for k = 1:2
     A(ell(j,k),ell(j,k)) = A(ell(j,k),ell(j,k))+1/h(j);
     integranda = @(xx) gfun(xx).*...
	 ((2-k)*(x(j+1)-xx)/h(j)+(k-1)*(xx-x(j))/h(j));
     gexact(ell(j,k)) = gexact(ell(j,k))+quadl(integranda,x(j),x(j+1),1e-12);
     g(ell(j,k)) = g(ell(j,k))+gell(j)*h(j)/2; % barycentric formula
     for i = k+1:2
       A(ell(j,k),ell(j,i)) = A(ell(j,k),ell(j,i))-1/h(j);
       A(ell(j,i),ell(j,k)) = A(ell(j,k),ell(j,i));      
     end
    end
  end
% Dirichlet b.c.
  A(1,1) = 1e16;
  gexact(1) = 0;
  g(1) = 0;
  A(m,m) = 1e16;
  gexact(m) = 0;
  g(m) = 0;
% solution
  u = A\gexact;
  ubar = A\g;
% linear interpolation
  uu = interp1(x,u,xx,'linear');
  uubar = interp1(x,ubar,xx,'linear');
% approximation of the L^2 norm (by trapezoidal rule on xx)
  errL2(counter) = norm(uu-uexact(xx))/sqrt(mm-1);
  errL2bar(counter) = norm(uubar-uexact(xx))/sqrt(mm-1);
%  errL2(counter) = norm(u-uexact(x))/sqrt(m-1); interpolant!
%  errL2bar(counter) = norm(ubar-uexact(x))/sqrt(m-1);
end
figure
loglog(mrange,errL2bar,'*g',mrange,mrange.^(-2),mrange,errL2,'ok')
legend('error in L2 norm','order 2','error in L2 norm (exact rhs)')
axis([10,1e3,1e-7,1e-1])
xlabel('m')
ylabel('error')
print('-depsc','fem1dnonsmooth.eps')