如何正确求解线性时变系统 附Matlab代码排查实现错误

核心错误点 #

• 状态转移矩阵括号错位：积分项中矩阵指数的参数存在语法级错误，原代码写为expm(A0*((k+1)*Ts)-u)，等价于对矩阵A0*(k+1)*Ts的所有元素减去标量u，完全违背状态转移矩阵expm(A0*((k+1)*Ts - u))的数学定义，是结果完全偏离预期的核心原因。
• 时变参数时间索引偏移：原代码直接用循环序号k作为时变参数g的时间自变量，未和实际仿真时间t对齐，g值全程存在固定时间偏移；同时循环多执行1次，计算区间超出预设的0~800仿真范围。
• 状态存储逻辑不匹配：原循环每次迭代先存储当前状态再推进计算，配合多余的循环步数，最终状态序列长度和时间轴t不匹配，绘图时会出现维度错位。
• 不必要的数值积分：单步内A、B取常数的零阶保持假设下，输入项积分存在解析解，用增广矩阵求矩阵指数的方式替代数值积分，精度更高、运算速度更快。

修正后可运行代码 #

r0 = 0.05;
L = 0.1;
d = 0.005;
w0 = 1.5;

Ts = 10;
t = 0:Ts:800;
N = length(t);
x0 = [0 0 1]';
x_value = zeros(3, N); % 预分配内存提升运行效率
x_value(:, 1) = x0; % 存入t=0时刻初始状态

for k = 1:N-1 % 共N-1个迭代步，从t(k)推进到t(k+1)
    t_cur = t(k);
    % 计算当前步起点的时变参数g，零阶保持假设下单步内g为常数
    g_cur = (2*r0*L*sinh((d*t_cur)/2)) / (d*cosh((d*t_cur)/2) + L*sinh((d*t_cur)/2));
    A0 = [-0.5*g_cur,  -w0,        0;
           w0,        -0.5*g_cur,  0;
           0,         0,          -g_cur];
    B0 = [0; 0; -g_cur];
    % 增广矩阵法一步求解状态，替代数值积分，无离散截断误差
    M = [A0, B0; zeros(1,3), 0];
    expM = expm(M*Ts);
    x_next = expM(1:3, 1:3)*x0 + expM(1:3, 4);
    % 存储状态并更新迭代初值
    x_value(:, k+1) = x_next;
    x0 = x_next;
end

plot(t, x_value, 'r-', 'linewidth', 1);
legend('x_1','x_2','x_3');
grid on;

若需要更高求解精度，可直接使用MATLAB内置的ode45等变系数微分方程求解器，避免大步长零阶保持离散引入的模型误差。

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