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')