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$. 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.csv');
disp(T1);
figure(1)
loglog(errorL2);
grid on;
grid minor;
end
% function to compute the L2 for an input of u, n and f source function
function errL2 = compute_error(u,n,f)
%solve the Poisson problem first.
[xgrid,h]=generate1d_uniform(n);
uh=fem1d(n,f);
% exact solution
uexact=u(xgrid);
% compute the error
errL2 = (h^1/2)*norm(uexact(:)-uh(:),2);
end