MATLAB roots函数实战指南：从系数构建到工程应用的全解析
引言
在工程计算和科学研究的道路上，多项式方程求解是一个无法绕开的基础课题。无论是控制系统的稳定性分析，还是信号处理的滤波器设计，甚至是金融模型的参数估计，都离不开多项式求根这一核心操作。MATLAB作为工程领域的标准工具，其内置的roots函数提供了一种高效可靠的解决方案。
但许多初学者在使用roots函数时常常陷入各种误区：系数顺序排列错误导致结果完全偏离预期，面对复数根时不知所措，或者对结果的准确性心存疑虑。这些问题不仅影响工作效率，更可能造成严重的工程判断失误。本文将从一个实践者的角度，通过大量工程案例和实用技巧，带你深入掌握roots函数的正确使用方式。
1. 理解多项式与系数向量的本质
1.1 多项式在MATLAB中的表示方法
在MATLAB中，多项式通过系数向量来表示，这是一个看似简单却容易出错的环节。关键在于理解系数排列的顺序规则：
- 降幂排列：系数向量必须按照变量的降幂顺序排列
- 完整包含所有幂次：即使某些幂次的系数为零，也必须保留其位置
- 向量方向无关：行向量或列向量均可接受
例如，多项式3x² + 2x - 2对应的系数向量为[3 2 -2]，而x⁴ - 1则表示为[1 0 0 0 -1]。
1.2 常见错误与验证方法
初学者常犯的错误包括：
- 系数顺序颠倒：将
[a b c]误写为[c b a] - 忽略零系数项：对于x³ + 2x - 1误写为
[1 2 -1]（正确应为[1 0 2 -1]） - 错误理解高次项：混淆x⁵和x⁴的系数位置
验证系数向量正确性的简单方法：
% 验证示例：检查多项式表示是否正确 p = [1 0 2 -1]; % 代表x³ + 2x - 1 x = 2; % 测试点 polyval(p, x) % 计算结果应为8 + 0 + 4 -1 = 112. roots函数的核心用法与深度解析
2.1 基础语法与参数说明
roots函数的基本调用形式极为简洁：
r = roots(p)其中：
p：多项式系数向量（如前所述规则）r：返回的根向量，总是列向量形式
2.2 实数根与复数根的处理
MATLAB不会区分实数根和复数根，所有根都以相同方式返回。对于实系数多项式，复数根总是以共轭对形式出现。
工程案例：RLC电路的特征方程求解
考虑一个RLC串联电路，其特征方程为： s² + (R/L)s + 1/(LC) = 0
% 参数设置 R = 100; % 电阻(Ohm) L = 0.1; % 电感(H) C = 1e-6; % 电容(F) % 构建系数向量 p = [1, R/L, 1/(L*C)]; % 求根 circuit_roots = roots(p) % 结果解释 if all(imag(circuit_roots) ~= 0) disp('电路处于欠阻尼状态，存在振荡') elseif any(circuit_roots > 0) disp('电路不稳定，存在发散风险') else disp('电路处于过阻尼状态，无振荡') end2.3 结果验证与误差分析
为确保roots结果的可靠性，可以采用以下验证方法：
- polyval验证：计算多项式在各根处的值，理论上应为零
- poly函数反向验证：使用
poly函数从根重建多项式系数 - 相对误差分析：比较原始系数与重建系数的差异
% 验证示例 p = [1 -3 3 -1]; % (x-1)^3 = x³ - 3x² + 3x -1 r = roots(p); reconstructed_p = poly(r); % 计算相对误差 relative_error = norm(p - reconstructed_p)/norm(p); disp(['相对误差：', num2str(relative_error)])3. 高阶应用技巧与工程实践
3.1 缺项多项式的处理方法
对于存在缺项的高阶多项式，必须严格保留零系数位置。例如，x⁵ + 2x² -1的系数向量应为[1 0 0 2 0 -1]。
机械振动分析案例：
考虑一个五自由度振动系统的特征方程： λ⁵ - 6λ⁴ + 10λ³ - 6λ² + λ = 0
% 注意λ的零次项系数为0 p = [1 -6 10 -6 1 0]; system_roots = roots(p); % 可视化根分布 figure plot(real(system_roots), imag(system_roots), 'ro', 'MarkerSize', 10) grid on title('系统特征根分布') xlabel('实部') ylabel('虚部')3.2 与poly函数的联合应用
poly函数是roots的逆运算，可用于：
- 从已知根重建多项式
- 验证
roots计算结果 - 构造特定根分布的系统
% 设计一个滤波器，要求零点在特定位置 desired_zeros = [-1, -0.5+0.5i, -0.5-0.5i]; filter_numerator = poly(desired_zeros); % 计算频率响应 freqz(filter_numerator, 1)3.3 数值稳定性与病态问题
对于高阶多项式（通常n>20），roots函数可能会遇到数值稳定性问题。此时可考虑：
- 平衡处理：对多项式进行变量替换或缩放
- 分段求解：将高次多项式因式分解为低次多项式乘积
- 替代算法：使用更稳定的求根算法如Jenkins-Traub
% 病态多项式示例：威尔金森多项式 p = poly(1:20); % 根为1到20的整数 r = roots(p); % 观察最大误差 max_error = max(abs(r - (1:20)')); disp(['最大根误差：', num2str(max_error)])4. 工程实战：控制系统稳定性分析
4.1 闭环系统特征方程求解
在控制系统中，稳定性由特征方程的根决定。左半平面的根代表稳定模态。
% 给定系统开环传递函数 num = [1]; den = [1 3 2 0]; % s³ + 3s² + 2s % 计算单位反馈下的闭环特征方程 K = 5; % 控制器增益 closed_loop_den = den + K*[zeros(1,length(den)-length(num)), num]; poles = roots(closed_loop_den); % 稳定性判断 if all(real(poles) < 0) disp('系统稳定') else disp('系统不稳定') unstable_poles = poles(real(poles) >= 0); disp(['不稳定极点：', num2str(unstable_poles')]) end4.2 根轨迹分析中的关键点计算
根轨迹与虚轴的交点可通过求解特定方程得到。
% 示例系统 num = [1]; den = conv([1 1], [1 2]); % 构造辅助方程求分离点 syms s aux_eqn = diff(poly2sym(num, s)*poly2sym(den, s)); break_points = double(roots(sym2poly(aux_eqn))); % 筛选实数分离点 valid_break_points = break_points(imag(break_points) == 0 & ... isreal(break_points)); disp(['实分离点：', num2str(valid_break_points')])4.3 多物理场耦合系统分析
对于复杂耦合系统，特征方程可能包含不同物理量纲的参数。
% 机电系统参数 J = 0.01; % 转动惯量(kg·m²) b = 0.1; % 阻尼系数(N·m·s) K = 0.5; % 弹簧常数(N·m/rad) R = 1; % 电阻(Ohm) L = 0.5; % 电感(H) Km = 0.3; % 电机常数(N·m/A) % 构建特征方程 p = [L*J, (R*J + L*b), (R*b + K*L + Km^2), K*R]; system_roots = roots(p); % 时域响应特性分析 damping_ratios = -real(system_roots)./abs(system_roots); natural_frequencies = abs(system_roots); disp(['阻尼比：', num2str(damping_ratios')]) disp(['自然频率(rad/s)：', num2str(natural_frequencies')])
