应用有限元法求解一组一阶常微分方程组的解的举例以及MATLAB代码实现

假设我们要求解以下的一组一阶常微分方程组：

$$
\begin{cases}
\frac{dx}{dt} = 2x + y \\
\frac{dy}{dt} = x + 3y
\end{cases}
$$

我们可以采用有限元法进行求解，具体步骤如下：

1. 将时间轴分为若干个小区间，每个小区间内的解使用一个多项式进行逼近。

2. 对于每个小区间 $[t_n, t_{n+1}]$，我们使用一个二次多项式进行逼近：

$$
x(t) = a_n + b_n(t-t_n) + c_n(t-t_n)^2 \\
y(t) = d_n + e_n(t-t_n) + f_n(t-t_n)^2
$$

其中 $a_n, b_n, c_n, d_n, e_n, f_n$ 都是待求解的系数。

3. 将上述多项式代入原方程组中，得到在每个小区间内的一组线性方程组。

4. 将相邻两个小区间内的解进行匹配，得到另外一组线性方程组。

5. 将上述两组线性方程组联立起来，得到总的线性方程组，可以使用 MATLAB 中的线性方程求解函数求解。

下面是 MATLAB 代码实现：

% 定义初始条件和时间范围
x0 = 1;
y0 = -1;
tspan = [0, 10];

% 定义有限元方法的参数
N = 100; % 将时间轴分为 N 个小区间
h = (tspan(2) - tspan(1)) / N; % 小区间的长度

% 定义多项式系数向量
a = zeros(N+1, 1);
b = zeros(N+1, 1);
c = zeros(N+1, 1);
d = zeros(N+1, 1);
e = zeros(N+1, 1);
f = zeros(N+1, 1);

% 初始化多项式系数向量
a(1) = x0;
d(1) = y0;

% 定义线性方程组的系数矩阵和右端向量
A = zeros(2*N, 2*N);
b = zeros(2*N, 1);

% 构造线性方程组
for n = 1:N
    % 在每个小区间内使用一个二次多项式进行逼近
    syms t;
    x = a(n) + b(n)*(t - n*h) + c(n)*(t - n*h)^2;
    y = d(n) + e(n)*(t - n*h) + f(n)*(t - n*h)^2;
    
    % 将多项式代入方程组中，得到在每个小区间内的一组线性方程组
    eq1 = diff(x) - 2*x - y;
    eq2 = diff(y) - x - 3*y;
    eqs = [eq1, eq2];
    vars = [b(n), c(n), e(n), f(n+1)];
    [A1, b1] = equationsToMatrix(eqs, vars);
    
    % 将每个小区间内的线性方程组加入总的线性方程组中
    A(n*2-1:n*2, n*2-1:n*2+2) = A1;
    b(n*2-1:n*2) = b1;
    
    % 匹配相邻两个小区间内的解，得到另外一组线性方程组
    eq1 = a(n+1) - (a(n) + b(n)*h + c(n)*h^2);
    eq2 = d(n+1) - (d(n) + e(n)*h + f(n)*h^2);
    A2 = [0, 0, 0, 1];
    b2 = [eq1; eq2];
    
    % 将相邻小区间内的线性方程组加入总的线性方程组中
    A(n*2+1:n*2+2, n*2-1:n*2+2) = A2;
    b(n*2+1:n*2+2) = b2;
end

% 求解线性方程组
x = A \ b;

% 提取出 x 和 y 的值
x_vals = x(1:2:end);
y_vals = x(2:2:end);

% 绘制 x 和 y 的变化曲线
plot(0:h:tspan(2), x_vals, 'r', 0:h:tspan(2), y_vals, 'b');
legend('x(t)', 'y(t)');
xlabel('t');
ylabel('x(t), y(t)');

运行以上 MATLAB 代码，即可得到该一组一阶常微分方程组的数值解。