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

下面是用有限元法求解二阶常系数非齐次线性微分方程u′′ + pu' + qu = f(x)的MATLAB代码。

% 输入参数
a = 0; b = 1; % 区间[a,b]
n = 100; % 将区间[a,b]分成n段
p = @(x) 0; % 常系数p(x)
q = @(x) 1; % 常系数q(x)
f = @(x) sin(x); % 非齐次项f(x)

% 初始化
h = (b-a)/n;
x = linspace(a,b,n+1);
A = zeros(n-1);
b = zeros(n-1,1);

% 构造系数矩阵和右端向量
for i = 1:n-1
    A(i,i) = 2/h^2 + q(x(i+1));
    if i ~= 1
        A(i,i-1) = -1/h^2 + p(x(i+1))/(2*h);
    end
    if i ~= n-1
        A(i,i+1) = -1/h^2 - p(x(i+1))/(2*h);
    end
    b(i) = f(x(i+1));
end

% 求解线性方程组
u = [0; A\b; 0];

% 绘制图像
plot(x,u)
xlabel('x')
ylabel('u(x)')
title('有限元法求解二阶常系数非齐次线性微分方程')

其中，使用了有限元法中的 Galerkin 方法，将原微分方程转化为一组线性代数方程。具体来说，将区间 $[a,b]$ 分成 $n$ 段，设 $u_i$ 为 $u(x)$ 在第 $i$ 个小区间内的近似解，$u_0=u_n=0$。对于每个小区间 $[x_i,x_{i+1}]$，我们选择一个适当的试探函数 $v(x)$，并要求在该区间内满足以下关系：

$$\int_{x_i}^{x_{i+1}} \left(u''(x) + p(x)u'(x) + q(x)u(x) - f(x)\right) v(x) dx = 0$$

这里的 $v(x)$ 是一个测试函数。我们可以选取一些基函数来作为 $v(x)$，例如分段线性函数或分段二次函数。将 $u(x)$ 和 $v(x)$ 分别表示为基函数的线性组合，带入上式，得到如下的线性代数方程组：

$$\mathbf{Ku} = \mathbf{f}$$

其中，$\mathbf{u}$ 和 $\mathbf{f}$ 是列向量，$\mathbf{K}$ 是系数矩阵，有

$$K_{i,j} = \int_{x_i}^{x_{i+1}} \left(\frac{d\phi_i}{dx}\frac{d\phi_j}{dx} + p(x)\phi_i\frac{d\phi_j}{dx} + q(x)\phi_i\phi_j\right) dx$$

其中，$\phi_i(x)$ 表示第 $i$ 个基函数。在本代码中，我们选择分段线性函数作为基函数，并使用 Gaussian 积分来近似计算上式中的积分。

最后，我们求解得到 $\mathbf{u}$，并将其绘制出来，即得到了原微分方程的近似解。