如何在Octave中绘制常微分方程组解x(1)与x(2)的曲线？

我已求解得到一个常微分方程组的解矩阵x（2×1符号矩阵），具体形式如下：

x = (sym 2×1 matrix)

  ⎡            -331⋅π⋅t                 -61⋅π⋅t  ⎤
  ⎢            ─────────                ──────── ⎥
  ⎢               2299               2   10835   ⎥
  ⎢ 689986800⋅ℯ            89076600⋅π ⋅ℯ         ⎥
  ⎢ ──────────────────── + ───────────────────── ⎥
  ⎢       1307333                 12173551       ⎥
  ⎢                                              ⎥
  ⎢             -331⋅π⋅t                 -61⋅π⋅t ⎥
  ⎢             ─────────                ────────⎥
  ⎢                2299                   10835  ⎥
  ⎢  206492400⋅ℯ            139429800⋅π⋅ℯ        ⎥
  ⎢- ──────────────────── + ─────────────────────⎥

以下是复现解x计算过程的代码：

pkg load symbolic;
A = [-0.4, 0.4; 0.05, -0.07];
[M, Lambda] = eig(A);
t = sym('t');
Phi = M * expm(Lambda*t)*inv(M);
x = Phi * [600; 0]

绘制x(1)和x(2)曲线的步骤 #

因为x是符号矩阵，需要先将其转换为可数值计算的函数，再生成时间点计算数值并绘图：

1. 将符号解转换为匿名函数

% 提取x的两个分量并转为匿名函数
x1_fun = matlabFunction(x(1));
x2_fun = matlabFunction(x(2));

2. 生成时间数组
   根据解的衰减特性，选择合适的时间范围（示例中取0到50）：

t_start = 0;
t_end = 50;
t_values = linspace(t_start, t_end, 1000); % 生成1000个均匀采样点

3. 计算各分量的数值

x1_values = x1_fun(t_values);
x2_values = x2_fun(t_values);

4. 绘制曲线

figure;
plot(t_values, x1_values, 'r-', 'LineWidth', 2, 'DisplayName', 'x(1)');
hold on;
plot(t_values, x2_values, 'b--', 'LineWidth', 2, 'DisplayName', 'x(2)');
xlabel('时间 t');
ylabel('解的分量');
title('常微分方程组解的分量随时间变化');
legend;
grid on;

注意事项 #

• 确保已通过pkg load symbolic加载符号计算包；
• 可根据解的实际衰减速度调整t_end的取值，让曲线变化趋势完整展示；
• 若需要更平滑的曲线，可增加linspace的采样点数。

内容的提问来源于stack exchange，提问作者Pavel Skipenes