0
$\begingroup$

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: enter image description here

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. enter image description here

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 
$\endgroup$
6
  • $\begingroup$ What are you plotting in the second graph? $\endgroup$ – David Apr 5 at 8:06
  • 3
    $\begingroup$ For this particular case of f(x)=-2 for which the solution is a quadratic function, Galerkin FEM gives nodally exact solutions with P1 (linear) elements. That is why your error is already at the machine precision and it does not change as you change the element size. If you want to see the second-order convergence, then you should use a more complex source term. $\endgroup$ – Chenna K Apr 5 at 8:26
  • $\begingroup$ @David I plot the log-scale error in the second plot. The x-axis should be the element size, but I don't plot details because I first wanted to see the order of convergence. $\endgroup$ – Dong Le Apr 5 at 15:53
  • 3
    $\begingroup$ To clarify @ChennaK's answer: what you are computing is the vector $\ell_2$ norm and you are comparing nodal values. It is not surprising that you get such small errors due to the interpolation properties of continuous P1 elements. If you compute the normal $L_2$ error using a quadrature rule, it should give you the correct rates. $\endgroup$ – Abdullah Ali Sivas Apr 5 at 20:30
  • 2
    $\begingroup$ You can take the average of i-th and (i+1)-st indices. That should work. $\endgroup$ – Abdullah Ali Sivas Apr 6 at 0:41

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.