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

下面是用 Matlab 实现有限元法求解二阶常系数非齐次线性微分方程的源代码，其中 P 和 Q 分别为常数系数，f(x)为非齐次项，a、b 为区间端点。

function [x,y] = FEM(P,Q,f,a,b,n)
% P,Q: 常系数
% f: 非齐次项
% a,b: 区间端点
% n: 划分的网格数
h = (b-a)/n;
% 生成节点及单元
x = linspace(a,b,n+1);
T = zeros(n,2);
for i = 1:n
    T(i,:) = [i,i+1];
end
% 组装刚度矩阵和载荷向量
A = zeros(n+1,n+1);
F = zeros(n+1,1);
for i = 1:n
    % 计算局部刚度矩阵和载荷向量
    k = [1/h,-1/h;-1/h,1/h];
    f1 = f(x(i));
    f2 = f(x(i+1));
    b = [f1*h/2+f2*h/6;f2*h/2+f1*h/6];
    % 组装到全局刚度矩阵和载荷向量
    A(T(i,:),T(i,:)) = A(T(i,:),T(i,:)) + P/h^2*k + Q/2/h*k;
    F(T(i,:)) = F(T(i,:)) + b;
end
% 处理边界条件
A(1,:) = 0;
A(1,1) = 1;
F(1) = 0;
A(n+1,:) = 0;
A(n+1,n+1) = 1;
F(n+1) = 0;
% 解方程
y = A\F;
end

使用示例：

P = 1; Q = 1; f = @(x)sin(x); a = 0; b = pi/2; n = 10;
[x,y] = FEM(P,Q,f,a,b,n);
plot(x,y);

其中，P=1，Q=1，f(x)=sin(x)，a=0，b=π/2，n=10，代表求解 y''+y'+y=sin(x)，区间为 [0,π/2]，将区间分成10个网格的解。