我写了一个 Matlab 代码来求解方程$-u'' = f$条件在域上。我用测试了代码。我用元素之一检查绘图，结果同意：$u(0) = u'(1) = 0$$x \in [0,1]$$f(x) = -2, \forall x\in [0,1]$$n=8$$u = u^{h}$

但是，我正在测试问题。我得到的错误图没有给我的结果，其中是固定步长。我正在尝试查找是否存在错误。为什么我的错误在和之间有很大的变化？这可能是我的代码中的一些小东西，但由于某种原因我无法弄清楚。提前致谢！！！！此外，我使用统一网格，因此是固定步长。 $L_{2}$$n=8,16,32,64$$O(\Delta x^{2})$$\Delta x$$n=16$$n=32$$[0,1]$$\Delta x (h)$

这是我用 Matlab 解决的代码和测试收敛性的驱动程序：

% 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