[Matlab科学计算] 有限元法求二阶常系数非齐次线性微分方程u′′ + pu' + qu = f(x)源代码，区间为a≤x≤b

以下是使用有限元法求解二阶常系数非齐次线性微分方程的Matlab程序。该程序使用线性三节点元素进行离散化，并使用高斯积分计算元素刚度矩阵和元素载荷向量。

% 输入参数
a = 0; % 区间左端点
b = 1; % 区间右端点
p = 1; % 二阶导数系数
q = 1; % 一阶导数系数
f = @(x) sin(x); % 非齐次项函数
n = 10; % 离散化节点数

% 初始化
h = (b-a)/(n-1); % 步长
A = zeros(n-2); % 系数矩阵
b = zeros(n-2,1); % 右端向量

% 计算元素刚度矩阵和载荷向量
K = [7/3 -4/3 1/3; -4/3 8/3 -4/3; 1/3 -4/3 7/3]; % 线性三节点元素刚度矩阵
F = [1/3; 4/3; 1/3]; % 线性三节点元素载荷向量

% 组装全局刚度矩阵和载荷向量
for i = 1:n-2
    % 计算元素左右节点位置
    x1 = a + (i-1)*h;
    x2 = a + i*h;
    x3 = a + (i+1)*h;
    
    % 计算元素刚度矩阵和载荷向量
    Ke = (p*h/30 + q*h^3/840)*K - 2*q*h^3/420*eye(3);
    Fe = h/6*F;
    Fe(1) = Fe(1) + f(x1)/2;
    Fe(2) = Fe(2) + f((x1+x2)/2);
    Fe(3) = Fe(3) + f(x2)/2;
    
    % 组装全局刚度矩阵和载荷向量
    A(i,i) = A(i,i) + Ke(1,1);
    A(i,i+1) = A(i,i+1) + Ke(1,2);
    A(i,i+2) = A(i,i+2) + Ke(1,3);
    if i<n-2
        A(i+1,i) = A(i+1,i) + Ke(2,1);
        A(i+1,i+1) = A(i+1,i+1) + Ke(2,2);
        A(i+1,i+2) = A(i+1,i+2) + Ke(2,3);
    end
    if i<n-3
        A(i+2,i) = A(i+2,i) + Ke(3,1);
        A(i+2,i+1) = A(i+2,i+1) + Ke(3,2);
        A(i+2,i+2) = A(i+2,i+2) + Ke(3,3);
    end
    b(i) = b(i) + Fe(1);
    if i<n-2
        b(i+1) = b(i+1) + Fe(2);
    end
    if i<n-3
        b(i+2) = b(i+2) + Fe(3);
    end
end

% 处理边界条件
A(1,:) = [];
A(:,1) = [];
b(1) = b(1) - (7/6*p/h + 1/6*q*h^2)*f(a);
A(end,:) = [];
A(:,end) = [];
b(end) = b(end) - (7/6*p/h + 1/6*q*h^2)*f(b);

% 解方程
u = [0; A\b; 0];

% 绘制解析解和数值解的比较图
x = linspace(a,b,100);
u_exact = -sin(x)/2;
plot(x,u_exact,'-',x,u,'o');
xlabel('x');
ylabel('u(x)');
legend('Exact solution','Numerical solution');

需要注意的是，该程序中使用的是线性三节点元素，因此需要至少有三个节点才能进行离散化。另外，需要特别处理边界条件，以满足方程的要求。