MATLAB多项式运算实战：从求值求根到数据拟合

1. 多项式运算：从基础概念到MATLAB实战

如果你正在学习工程、物理、数学或者任何涉及数据分析的领域，那么“多项式”这个概念你一定绕不开。它看起来简单——不就是像x^2 + 2x + 1这样的表达式吗？但在实际应用中，无论是信号处理中的滤波器设计，控制系统中的传递函数分析，还是机器学习中的特征拟合，多项式都扮演着核心角色。今天，我们不谈枯燥的教科书理论，而是从一个实践者的角度，聊聊如何在MATLAB这个强大的工具里，高效、准确地“玩转”多项式。你会发现，从简单的求值、求导，到复杂的求根和拟合，MATLAB提供了一套极其顺手的工具箱，能让你把抽象的数学公式瞬间变成可视化的结果和可验证的代码。

2. 多项式在MATLAB中的核心表达与基础操作

在动手之前，我们必须统一“语言”。在MATLAB的世界里，多项式有它自己的一套“语法”，理解这套语法是高效操作的前提。

2.1 多项式的向量表示法

这是MATLAB处理多项式最核心、也最需要习惯的一点。MATLAB不直接存储x^3 - 2x - 5这样的符号表达式（除非使用符号计算工具箱）。相反，它用一个行向量来存储多项式的系数，并且顺序是固定的：从最高次幂到常数项。

举个例子，对于多项式：P(x) = 4x^5 + 3x^4 - 2x^2 + x + 7

它在MATLAB中应该表示为：p = [4, 3, 0, -2, 1, 7]

注意看，x^3的系数是0，这个位置绝对不能省略！你必须用0来占位，以确保向量中系数的索引位置与x的幂次严格对应。向量p的第一个元素4对应x^5，最后一个元素7对应常数项。

实操心得：新手最容易犯的错误就是漏写零系数。一个快速检查的方法是，你的系数向量长度应该等于（多项式最高次数 + 1）。上面例子中，最高次是5，向量长度就是6。养成这个习惯，能避免很多诡异的计算错误。

2.2 基础三板斧：求值、求导与四则运算

一旦用向量定义好了多项式，一系列基础操作就变得异常简单。

1. 多项式求值：polyval函数这是最常用的函数之一。给定一个多项式系数向量p和一组自变量x，polyval(p, x)会返回对应的函数值y。

p = [1, -5, 6]; % 代表 x^2 - 5x + 6 x_points = [0, 2, 3, 4]; y_values = polyval(p, x_points); disp(y_values); % 输出：6 0 0 2

这行代码计算了多项式在 x=0,2,3,4 四个点的值。在需要绘制多项式函数图像时，polyval是必不可少的。

2. 多项式求导：polyder函数这就是热搜词里的POLYDER。它直接对系数向量进行操作，返回求导后新多项式的系数向量。

p = [2, 0, -4, 3]; % 代表 2x^3 - 4x + 3 dp = polyder(p); % 求导 disp(dp); % 输出：6 0 -4， 代表导数 6x^2 - 4

polyder还可以用于求两个多项式乘积或商的导数，语法分别为polyder(a,b)和[num, den] = polyder(b,a)（返回商函数的分子分母导数系数）。这在控制系统分析中求取传递函数的导数时非常有用。

3. 多项式乘法与除法：conv和deconv函数多项式的乘法对应于其系数向量的卷积运算。

a = [1, 2]; % x+2 b = [1, -3]; % x-3 c = conv(a, b); % 计算 (x+2)(x-3) disp(c); % 输出：1 -1 -6， 代表 x^2 - x - 6

而除法是乘法的逆运算，使用deconv，它返回商和余数。

[quotient, remainder] = deconv(c, a); % 用 c 除以 a disp(quotient); % 输出：1 -3， 即 b disp(remainder); % 输出：0， 能整除

3. 核心进阶：求根与由根构建多项式

多项式的“根”（或“零点”）是其灵魂所在。在工程上，根的分布直接决定了系统的稳定性、滤波器的频率特性等关键性质。

3.1 求取多项式根：roots函数

roots函数是另一个明星函数，也是热搜词ROOTS所指。它接收一个多项式系数向量p，返回一个包含其所有根（可能是实数或复数）的列向量。

p = [1, -5, 6]; % x^2 - 5x + 6 r = roots(p); disp(r); % 输出：3.0000 和 2.0000

对于高次多项式，roots函数通过计算伴随矩阵的特征值来求解，这是一种数值方法。这意味着对于非常高次或病态的多项式，结果可能存在数值误差。

注意事项：roots函数的结果是数值解。对于5次及以上的多项式，没有通用的求根公式，数值方法是唯一选择。此外，如果系数向量p的第一个元素是0，MATLAB会自动忽略开头的零，但中间的零系数必须保留。

3.2 由根构建多项式：poly函数

这是roots的逆过程。给定一个包含根的向量r，poly(r)会生成一个以这些根为零点的首一多项式（最高次项系数为1）的系数向量。

r = [2, -1, 3+4i, 3-4i]; % 设定一些根，包括复共轭根 p_reconstructed = poly(r); disp(p_reconstructed); % 输出一个系数向量，代表多项式 (x-2)(x+1)(x-(3+4i))(x-(3-4i)) % 展开后是一个实系数四次多项式。

这个功能极其强大。例如，在滤波器设计中，我们首先在复平面上确定所需的极点（根）位置，然后直接用poly函数得到系统函数的分子或分母多项式系数。

3.3 一个综合案例：验证根与系数的关系

让我们把polyval,roots,poly串起来，完成一个自我验证的闭环。

% 1. 原始多项式 p_original = [1, 0, -13, 0, 36]; % x^4 - 13x^2 + 36 % 2. 求根 r_calculated = roots(p_original); % 3. 由根重构多项式 p_reconstructed = poly(r_calculated); % 4. 比较（由于数值计算误差，直接比较可能不相等） % 我们可以比较在多个采样点上的值是否接近 x_test = linspace(-5, 5, 100); y_original = polyval(p_original, x_test); y_reconstructed = polyval(p_reconstructed, x_test); error = max(abs(y_original - y_reconstructed)); fprintf(‘最大误差为： %e\n‘， error); % 通常是一个非常小的数，如 1e-12 量级

这个案例清晰地展示了多项式在系数域和根域之间的转换，这是许多工程应用的数学基础。

4. 多项式拟合：从数据中寻找规律

前面我们讨论的是“已知多项式，分析其性质”。但现实中更常见的问题是：我有一堆离散的数据点(x, y)，我想找到一个多项式来近似描述它们之间的关系。这就是多项式拟合。

4.1 使用polyfit进行最小二乘拟合

MATLAB 提供了强大的polyfit函数。它的基本语法是p = polyfit(x, y, n)，其中n是你希望拟合的多项式的次数。函数会返回一个系数向量p，使得多项式P(x)在最小二乘意义下最优地逼近数据。

% 生成一些带有噪声的测试数据（假设背后是二次关系） x_data = linspace(0, 10, 30); y_true = 0.5 * x_data.^2 - 2 * x_data + 1; noise = randn(size(x_data)) * 2; % 加入高斯噪声 y_data = y_true + noise; % 进行二次多项式拟合 degree = 2; p_fit = polyfit(x_data, y_data, degree); % 绘制结果 x_fine = linspace(min(x_data), max(x_data), 200); y_fit = polyval(p_fit, x_fine); figure; scatter(x_data, y_data, ‘b‘， ‘DisplayName‘， ‘原始数据（含噪声）‘); hold on; plot(x_fine, y_fit, ‘r-‘， ‘LineWidth‘， 2, ‘DisplayName‘， sprintf(‘%d次拟合曲线‘， degree)); plot(x_data, y_true, ‘g--‘， ‘LineWidth‘， 1.5, ‘DisplayName‘， ‘真实函数‘); legend(‘Location‘， ‘best‘); xlabel(‘x‘); ylabel(‘y‘); title(‘多项式拟合示例‘); grid on;

运行这段代码，你会看到红色的拟合曲线如何穿过蓝色的噪声数据点，并逼近绿色的真实函数。

4.2 拟合阶数的选择：偏差与方差的权衡

选择拟合阶数n是一门艺术，也是实践中最容易出错的地方。

• 欠拟合 (Underfitting)：当n太小（例如用直线去拟合明显弯曲的数据），模型过于简单，无法捕捉数据中的潜在规律，导致偏差很大。在图中表现为拟合曲线远离大部分数据点。
• 过拟合 (Overfitting)：当n太大（例如用10次多项式拟合20个数据点），模型过于复杂，它不仅拟合了潜在规律，还“拟合”了噪声。这会导致方差很大，模型在训练数据上表现极好，但对新数据的预测能力（泛化能力）很差。在图中表现为拟合曲线剧烈波动，穿过每一个数据点。

实操心得：如何选择n？没有绝对答案，但可以遵循以下步骤：

1. 可视化：始终绘制拟合曲线与原始数据的对比图。这是最直观的方法。
2. 交叉验证：将数据分为训练集和测试集。用训练集拟合多项式，然后在测试集上计算误差（如均方误差MSE）。选择那个在测试集上误差最小的n。
3. 观察拟合系数：如果高阶项的系数绝对值非常小（例如小于1e-6），可能意味着该高阶项贡献不大。
4. 奥卡姆剃刀原则：在拟合效果相近时，优先选择阶数更低的简单模型。

4.3 评估拟合质量：polyval与残差分析

拟合完成后，我们需要量化评估其好坏。除了看图，计算残差是关键。

% 接上一段代码 % 计算在原始数据点上的拟合值 y_data_fit = polyval(p_fit, x_data); % 计算残差（误差） residuals = y_data - y_data_fit; % 绘制残差图 figure; subplot(2,1,1); plot(x_data, residuals, ‘ko-‘); xlabel(‘x‘); ylabel(‘残差‘); title(‘残差图‘); grid on; hline = refline(0,0); % 添加y=0参考线 hline.Color = ‘r‘; hline.LineStyle = ‘--‘; subplot(2,1,2); histogram(residuals, 20); xlabel(‘残差值‘); ylabel(‘频数‘); title(‘残差分布直方图‘); grid on;

一个“好”的拟合，其残差图应该呈现出随机散布的模式，没有明显的趋势或结构（如弯曲、漏斗形）。残差直方图应大致接近正态分布。如果残差图显示出明显的规律，说明当前的多项式模型可能不足以描述数据，需要考虑更复杂的模型或检查数据本身的问题。

5. 实战演练与疑难排坑指南

理论说得再多，不如动手踩一遍坑来得实在。下面我们通过一个综合性的实战案例，串联起大部分操作，并附上我多年来总结的“避坑指南”。

5.1 综合案例：设计一个简单低通滤波器的频率响应

在信号处理中，一个简单无限脉冲响应滤波器的系统函数常表示为多项式之比H(z) = B(z)/A(z)。其中A(z)的根（极点）决定了滤波器的特性。让我们设计一个截止频率相关的低通滤波器。

% 案例：设计一个二阶巴特沃斯低通滤波器（模拟原型） % 其极点位于左半平面单位圆上，角度与截止频率相关 wc = 0.5*pi; % 归一化截止频率 (0.5pi rad/sample) % 计算极点位置（二阶情况） theta = pi/4; % 对于二阶巴特沃斯，极点角度为45度 p1 = exp(1i * (pi + theta)); % 单位圆上的极点 p2 = exp(1i * (pi - theta)); % 注意：实际滤波器设计请使用 `butter`, `cheby1` 等专业函数，此处为演示多项式操作 % 由极点（根）构建分母多项式 A(z)（首一多项式） A = poly([p1, p2]); % A(z) = (z - p1)(z - p2) % 分子多项式 B(z) 通常为常数，使得直流增益为1 B = sum(A); % 一种简单设置，使 H(z=1)=1 % 计算频率响应 N_fft = 512; [H_freq, w] = freqz(B, A, N_fft); % w 是归一化角频率向量 % 绘制幅频响应 figure; plot(w/pi, 20*log10(abs(H_freq)), ‘b-‘， ‘LineWidth‘， 2); xlabel(‘归一化频率 (×π rad/sample)‘); ylabel(‘幅度响应 (dB)‘); title(‘二阶滤波器幅频响应（多项式构建演示）‘); grid on; xlim([0, 1]);

这个案例展示了如何从设定的根（极点）出发，利用poly函数构建系统函数的分母，进而分析系统特性。虽然实际工程中我们会用更专业的工具箱函数，但理解其背后的多项式操作至关重要。

5.2 常见问题与排查技巧实录

在长期使用中，我积累了一些典型问题的解决方法：

问题1：使用roots求高次多项式根时，结果出现微小虚部。

• 现象：对于一个实系数多项式，理论上非实根应以共轭对出现。但roots的数值计算可能使本应为实数的根产生一个非常小的虚部（如1.0e-14i）。
• 解决：使用real函数取实部。但需谨慎，应先判断虚部是否真的可忽略。

  r = roots(high_order_poly); % 方法：设定一个容差阈值 tolerance = 1e-10; r_real = real(r(abs(imag(r)) < tolerance)); disp(‘处理后的实根：‘); disp(r_real);

问题2：polyfit拟合高阶多项式时出现警告或结果异常（NaN）。

• 原因：这通常是因为病态范德蒙矩阵问题。当数据点x的跨度很大或拟合阶数n很高时，用于求解最小二乘问题的方程组条件数极大，数值不稳定。
• 解决：
  1. 中心化与缩放：将x数据减去其均值，并除以其标准差，再进行拟合。拟合后，需要对得到的多项式系数进行相应的逆变换。
  2. 使用正交多项式：更稳健的方法是使用polyfit的另一种形式（返回用于评估的结构体），或手动使用qr分解。
  3. 根本方法：重新考虑是否真的需要如此高的阶数。尝试降低拟合阶数n。

问题3：conv进行多项式乘法后，系数向量长度激增，导致后续polyval计算缓慢或不精确。

• 分析：两个次数分别为m和n的多项式相乘，结果次数为m+n。如果连续进行多次乘法，系数向量长度会呈指数增长，不仅计算慢，在x很大或很小时还容易导致浮点数上溢或下溢。
• 优化策略：
  • 如果最终目的是求值，考虑使用嵌套的polyval。例如，求P(x)*Q(x)在x处的值，优先计算polyval(P, x) .* polyval(Q, x)，而不是先conv再polyval。
  • 对于特别高次的多项式，考虑使用符号计算工具箱进行化简，或转换为其他表示形式（如切比雪夫多项式基）。

问题4：如何判断多项式是否有重根？

• 方法：重根在数值上表现为roots返回非常接近的值。但更严谨的方法是检查多项式与其导数的结式或求最大公因式。一个实用技巧是：

  p = [1, -3, 3, -1]; % (x-1)^3 r = roots(p); % 查看根之间的最小距离 min_distance = min(pdist(r)); if min_distance < 1e-8 warning(‘可能存在非常接近的根，即重根或临近根。‘); end % 精确方法：计算多项式与其导数的最大公因式（需要符号工具箱） % syms x; p_sym = poly2sym(p, x); dp_sym = diff(p_sym, x); % gcd_result = gcd(p_sym, dp_sym);

问题5：polyder求导后，如何求导函数在某点的值？

• 误区：不要试图从求导后的系数向量dp去反推导函数表达式再求值。
• 正确做法：有两种等效且更数值稳定的方法：
  1. 直接使用polyval对导数系数向量求值：derivative_at_x0 = polyval(dp, x0)。
  2. 使用多项式求值的底层原理——秦九韶算法（Horner‘s method），MATLAB的polyval已优化。对于高阶多项式，手动编写秦九韶算法求导数值可能更快，但polyval在大多数情况下已足够优秀。

处理多项式，尤其是用MATLAB这样的数值计算工具，核心在于理解其“系数向量”的表示本质，并时刻警惕数值计算中潜在的精度问题和稳定性问题。从简单的polyval、roots，到综合应用的polyfit，每一个函数都是将数学思想转化为工程实践的神奇桥梁。多动手尝试，多观察结果，特别是图形化结果，你会对多项式这个强大的数学工具有更深刻的直觉。最后，记住那句老话：模型越简单越好，在能用一次线性模型的时候，就不要为了追求“完美拟合”而强行使用十次多项式。