基于MATLAB ODE45求解结果的非线性系统特征值求解问询

关于直接代入X_answer到K函数句柄的问题 #

不能直接将t_ode×resSize的X_answer矩阵整体代入K(X)函数句柄，因为K的输入是单个时刻的位移向量（维度resSize×1），而非多时刻的位移矩阵。非线性系统的刚度K(X)随位移（时间）变化，其特征值是时变的——不同时刻的系统线性化特性不同，特征值也不同。

你需要针对单个时刻的位移向量计算K矩阵，再调用polyeig求解该时刻的特征值：

% 示例：遍历所有时刻，计算每个时刻的特征值
n_time = size(t_ode, 1);
lambda_all = cell(n_time, 1); % 存储每个时刻的特征值

for i = 1:n_time
    X_t = X_answer(i, :)'; % 取出第i时刻的位移向量（转置为列向量）
    K_t = K(X_t); % 计算该时刻的刚度矩阵
    [vector_t, lambda_t, cond_t] = polyeig(K_t, C, M);
    lambda_all{i} = lambda_t;
end

% 若只需分析关键时刻（如稳态时刻），直接取对应行即可
X_steady = X_answer(end, :)'; % 假设最后时刻为稳态
K_steady = K(X_steady);
[vector_steady, lambda_steady, cond_steady] = polyeig(K_steady, C, M);

基于雅可比矩阵的局部特征值求解 #

对于非线性二阶ODE：
$$M\ddot{X} + C\dot{X} + K(X) = F(t)$$
将其转化为一阶状态空间形式更便于计算雅可比矩阵：令状态向量$Y = \begin{bmatrix}X \ \dot{X}\end{bmatrix}$，则状态方程为：
$$\dot{Y} = \begin{bmatrix}\dot{X} \ M^{-1}(F(t) - C\dot{X} - K(X))\end{bmatrix}$$

雅可比矩阵$J$是$\dot{Y}$对$Y$的偏导矩阵（维度$2\text{resSize}×2\text{resSize}$），分块形式为：
$$J = \begin{bmatrix}0 & I \ -M^{-1}\frac{\partial K(X)}{\partial X} & -M^{-1}C\end{bmatrix}$$
其中：

• $0$是$\text{resSize}×\text{resSize}$零矩阵，$I$是$\text{resSize}×\text{resSize}$单位矩阵
• $\frac{\partial K(X)}{\partial X}$是刚度矩阵对位移向量$X$的偏导（若$K(X)$是显式函数，可手动推导或用MATLAB的jacobian函数符号计算）

计算出某时刻的雅可比矩阵$J$后，直接调用eig(J)即可得到该时刻系统的局部线性化特征值：

% 示例：计算某时刻的雅可比矩阵及特征值
X_t = X_answer(i, :)'; % 第i时刻位移
Xd_t = gradient(X_answer(:, i), t_ode); % 近似计算该时刻速度
Xd_t = Xd_t(i);

% 符号计算刚度矩阵的偏导（假设K(X)是符号表达式）
syms X_sym(resSize, 1);
K_sym = K(X_sym);
dKdX_sym = jacobian(K_sym, X_sym);
dKdX_t = double(subs(dKdX_sym, X_sym, X_t)); % 代入当前位移得到数值矩阵

% 构建雅可比矩阵
I_mat = eye(resSize);
O_mat = zeros(resSize);
M_inv = inv(M);
J = [O_mat, I_mat; -M_inv*dKdX_t, -M_inv*C];

% 计算特征值
lambda_jac = eig(J);

关于传递函数的说明 #

非线性系统不存在全局的传递函数，只有局部线性化后的传递函数——基于上述雅可比矩阵得到的线性状态空间模型，可转化为传递函数，其极点即为对应时刻的特征值。若需分析局部传递函数，可使用MATLAB的ss和tf函数：

% 基于雅可比矩阵构建线性状态空间模型
sys_ss = ss(J, zeros(2*resSize,1), eye(2*resSize), zeros(2*resSize,1));
% 转化为传递函数（以第一个状态输出为例）
sys_tf = tf(sys_ss(1, :));
% 提取极点（即特征值）
poles = pole(sys_tf);

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