用MATLAB求解微分方程组

MATLAB 可以方便地求解微分方程组，下面我们以一个二阶线性微分方程组为例进行说明。

假设有一个二阶线性微分方程组：

y1'' + 2y1' - 3y2 = 0

y2'' - 2y1' - 4y2 = cos(t)

其中，y1 和 y2 是待求解的函数，t 是自变量，cos(t) 是已知函数。这个方程组可以表示成向量形式：

Y' = F(t, Y)

其中，Y = [y1; y2]，F(t, Y) = [y1'; y2'] = [y1'; y2' - 2y1' - 4y2]，即：

Y' = [y1'; y2'] = [y1'; y2' - 2y1' - 4y2] = F(t, Y)

我们可以使用 MATLAB 中的 ode45 函数来求解这个方程组。

首先，我们需要将方程组表示成 MATLAB 可以识别的函数形式。定义一个函数文件（例如 myfun.m），其中包含方程组的函数表达式：

function dydt = myfun(t, y)
    dydt = zeros(2, 1);
    dydt(1) = y(2);
    dydt(2) = 2*y(2) + 3*y(1);
end

然后，我们可以调用 ode45 函数来求解方程组。调用格式为：

[t, Y] = ode45(@myfun, tspan, Y0)

其中，@myfun 是函数句柄，表示要求解的微分方程组对应的函数，tspan 是一个包含起始时间和终止时间的向量，Y0 是一个包含 y1 和 y2 的初始值的向量。

例如，我们可以设定起始时间为 0，终止时间为 10，初始值为 y1(0) = 0，y2(0) = 1，代码如下：

tspan = [0, 10];
Y0 = [0; 1];
[t, Y] = ode45(@myfun, tspan, Y0);

这样，t 就是一个包含时间节点的向量，Y 就是相应时间节点的 y1 和 y2 值。

如果方程组右侧还有一个非齐次项，我们只需要修改一下 myfun 函数即可。例如，如果右侧的非齐次项是 cos(t)，我们可以将 myfun 修改为：

function dydt = myfun(t, y)
    dydt = zeros(2, 1);
    dydt(1) = y(2);
    dydt(2) = 2*y(2) + 3*y(1) + cos(t);
end

然后再调用 ode45 函数即可。

需要注意的是，对于高阶微分方程组，我们需要将其转化为一阶微分方程组的形式，才能使用 ode45 函数求解。