Error $L_{2}$ convergence in Finite Element for Poisson Equation
I have written a Matlab code to solve the equation $-u'' = f$ with conditions $u(0) = u'(1) = 0$ on the domain $x \in [0,1]$. I tested the code with $f(x) = -2, \forall x\in [0,1]$. I check the plot with one of element $n=8$, and the result $u = u^{h}$ agrees:
However, I am testing the convergence error $L_{2}$ for the problem $n=8,16,32,64$. 935
The error plot I get is not giving me the result of $O(\Delta x^{2})$ where $\Delta x$ is the fixed step size. I am trying to find if there is a bug. Why is a big change in my error between $n=16$ and $n=32$? It is probably something small in my code, but for some reason I can't figure it out. Thanks in advance!!!! Furthermore, I use a uniform grid $[0,1]$ so $\Delta x (h)$ is a fixed step size.
This is the code that I solve with Matlab and a driver to test the convergence:
``` % script to generate the uniform mesh for Finite element: % dom = [a b] where a < b is the domain of our problem function [xgrid, h] = generate1d_uniform(n) h=1/(n-1); xgrid=linspace(0,1,n); end
% implement the solver for Finite Element for Poisson equation function uh = fem1d(n,f)
% generate the uniform mesh: [x,h] = generate1d_uniform(n); % load vector b = loadVector1D(x,h,f); A = StiffMat1D(x,h);
% since u0 is 0 Ae=A(2:end,2:end); fe=b(2:end); uhe=Ae\fe;
% adjust the boundary condition u(0)=0 uh=[0;uhe]; end
% function to assembly load vector b: % f is the source function; x is from the domain generate by 1d mesh function b = loadVector1D(x,h,f) n = length(x)-1; b = zeros(n+1,1);
for i = 1:n b(i) = b(i) + f(x(i))(h/2); b(i+1) = b(i+1) + f(x(i+1))(h/2); end
end
% function to load the stiffness matrix for A: % the local stiffness matrix is of the form Ak = 1/h [1 -1; -1 1] % adjust the A(1,1) = 1 and A(n+1,n+1) = 1 for boundary conditions in this % case u'(1) = 0 and u(0) = 0 function A = StiffMat1D(x,h) n = length(x)-1; Ak = spdiags(ones(n+1,1)[-1 2 -1],-1:1,n+1,n+1); % Adjust the boundary condition u'(1) = 0 Ak(1,1) = 1; Ak(end,end) = 1;
A = (1/h)Ak; end
%%% DRIVER for convergence error %%% function driver1 clear all; clc; format short; %number of elements set: nset=[8,16,32,64];
% source function in part (a) f = @(x) -2.*(0<=x & x<=1);
% exact solution: u = @(x) x.^2 - 2*x; %define the error L2 set to store the values: errorL2=zeros(size(nset)); hset=zeros(size(nset)); % iteration: niter=length(nset);
for i=1:niter
n=nset(i);
disp(n)
errorL2(i)=compute_error(u,n,f);
end
% write the table for error L2: T1=table(nset', errorL2'); T1.Properties.VariableNames ={'n elements', 'Error L2'}; writetable(T1,'TableError. 935