Graeffe根平方法：从原理到MATLAB实现，解决多项式求根数值难题
1. 从“暴力”到“优雅”：为什么我们需要Graeffe根平方法？
在工程计算、信号处理或者控制系统设计的日常里，多项式求根是个绕不开的活儿。你可能在MATLAB里随手敲个roots([1, -5, 6])，眨眼间就得到了方程x^2 - 5x + 6 = 0的根2和3。这太方便了，以至于我们很少去关心roots函数背后那台精密的“机器”是如何运转的。但当你面对一个次数高达几十甚至上百、系数跨度极大的多项式时，事情就开始变得棘手了。你可能会发现roots函数返回的结果里出现了诡异的复数，或者数值稳定性极差，一个微小的系数扰动就让根的位置“面目全非”。这时候，你需要的可能不是更强大的计算机，而是一种更稳健的算法。这就是我们今天要深入探讨的Graeffe根平方法。
Graeffe方法诞生于那个没有电子计算机的年代，它的核心思想充满了数学的巧思：与其直接和多项式的根“硬碰硬”，不如通过一系列巧妙的变换，把根的“信息”放大、分离，最终以一种非常数值稳定的方式把它们计算出来。它特别擅长处理那些根在复平面上模值相差很大的情况——而这恰恰是很多工程问题（比如滤波器设计、系统稳定性分析）中多项式根分布的典型特征。虽然如今它不再是MATLABroots函数的默认选择（后者基于更现代的“伴侣矩阵”特征值算法），但理解Graeffe方法，不仅能让你在遇到数值困难时多一个备选工具，更能深刻体会到数值分析中“问题转化”和“稳定性设计”的精妙之处。这篇文章，我将带你从零开始，手把手实现Graeffe方法，并深入剖析其原理、优势、局限以及在MATLAB环境下的实战技巧与避坑指南。
2. Graeffe方法的核心原理：把根的“距离”变成“尺度”
要理解Graeffe方法为什么有效，我们得先暂时忘掉那些复杂的迭代公式，从它的几何直觉入手。考虑一个简单的二次多项式P(x) = (x - r1)(x - r2) = x^2 + a1*x + a0。它的两个根是r1和r2。现在，我们构造一个新的多项式Q(x) = P(√x)*P(-√x)。这个操作看起来有点奇怪，但如果你动手算一下（或者信任我的推导），会发现一个惊人的事实：Q(x)的根，变成了r1^2和r2^2。
为什么这个事实如此重要？想象一下，如果|r1|和|r2|原本比较接近，比如一个是5，一个是6，那么直接区分它们对数值误差很敏感。但如果我们平方一下，它们就变成了25和36，差距从1拉大到了11！这个“差距”被显著放大了。Graeffe方法的精髓，就是反复执行这种“平方根”变换（实际上通过多项式系数操作实现，无需真的开方），在每一步迭代中，将原多项式根的平方（乃至更高次幂）作为新多项式的根。经过k次迭代后，新多项式的根就是原根的2^k次幂。
2.1 从原理到系数递推公式
当然，我们不可能每次都去显式地构造P(√x)*P(-√x)再展开。Graeffe的智慧在于，他找到了一种直接通过原多项式系数进行迭代，得到新多项式系数的递推关系。对于一个n次多项式：P(x) = x^n + a_{n-1}x^{n-1} + ... + a_1x + a_0
定义第m次迭代后的多项式系数为a_k^(m)（其中a_n^(m) ≡ 1）。那么Graeffe迭代的递推公式为：a_k^(m+1) = (-1)^k * [ (a_k^(m))^2 + 2 * Σ_{j=1}^{min(k, n-k)} (-1)^j * a_{k-j}^(m) * a_{k+j}^(m) ]
这个公式初看复杂，但其核心是交叉项的加减组合。它保证了迭代后的多项式P_m(x)的根恰好是原多项式根的2^m次幂。在实际编程中，我们正是利用这个公式，从初始系数a_k^(0)出发，一层层计算出a_k^(1),a_k^(2), ...。
注意：这个公式在系数为实数时，能有效处理实根和共轭复根。但对于更一般的复系数多项式，符号处理需要格外小心，通常建议先将算法在实系数情况下跑通。
2.2 迭代的终点：如何从放大后的根反推原根？
经过足够多次（比如s次）迭代后，多项式P_s(x)的根R_i已经变成了r_i^(2^s)。由于幂次极高，这些R_i的模长会呈现出极其夸张的差异（只要原根模长不完全相等）。此时，一个关键性质登场：对于模长占绝对优势的根，其对应的高次幂会“主导”其所在多项式的系数。
具体来说，考虑P_s(x)的常数项a_0^(s)。根据韦达定理，它等于所有根R_i的乘积（乘以(-1)^n）。如果有一个根R_1的模长远远大于其他所有根，那么|a_0^(s)| ≈ |R_1|。更进一步，如果前k个最大模长的根彼此之间也拉开了指数级差距，那么我们可以通过一系列系数比来逐级恢复它们。
最常用的恢复公式是基于相邻两次迭代的系数比。对于模长最大的根r_1，有：|r_1| ≈ lim_{s->∞} |a_{n-1}^{(s)} / a_{n-1}^{(s-1)}|^(1/2^s)而在有限次迭代后，一个更实用的近似是：|r_1|^(2^s) ≈ |a_{n-1}^{(s)} / a_{n-1}^{(s-1)}|实际上，我们通常对每个根r_i使用：|r_i| ≈ ( |a_{n-i}^{(s)} / a_{n-i}^{(s-1)}| )^(1/2^s)
得到模长后，还需要确定辐角（即根是实数还是复数，复数的角度是多少）。对于实系数多项式，实根的辐角是0或π（对应正负号），共轭复根的辐角互为相反数。辐角可以通过类似的思想，利用系数比值的符号和更精细的公式来恢复，但这一步的数值稳定性比模长计算要差，是实践中容易出问题的地方。
3. MATLAB实战：手把手实现一个稳健的Graeffe求根器
理解了原理，我们就在MATLAB里把它实现出来。我们的目标是写一个函数graeffeRoots(p, s)，其中p是多项式系数向量（按降幂排列），s是指定的迭代次数。
3.1 基础迭代框架与系数计算
首先，我们处理输入并初始化。一个健壮的程序应该能处理用户输入的行向量或列向量，并自动去除最高次项系数可能存在的零（虽然理论上应该是1）。
function roots_approx = graeffeRoots(p, max_iter) % 输入: p - 多项式系数向量，例如 [1, -5, 6] 代表 x^2 - 5x + 6 % max_iter - Graeffe迭代次数，通常8-12次足够 % 输出: roots_approx - 计算得到的根（复数向量） % 1. 预处理系数：确保是行向量，并归一化（使最高次项系数为1） p = p(:).'; % 转为行向量 if p(1) ~= 1 p = p / p(1); % 归一化 end n = length(p) - 1; % 多项式次数 % 2. 初始化系数矩阵，每一行存储一次迭代的系数 % a(k, i) 对应第k次迭代中 x^(n-i) 的系数？不，更直观的：a_coeffs(iter, idx) % 我们让 a_coeffs(iter, 1) 对应常数项 a0， a_coeffs(iter, end) 对应 x^n 系数（恒为1） % 但Graeffe递推公式通常索引从0到n。在MATLAB中，我们用索引 1 代表 a0。 a_prev = p(end:-1:1); % 初始系数 [a0, a1, ..., a_{n-1}, 1]？ 注意我们的p是降幂。 % 实际上，p = [1, a_{n-1}, ..., a0]，所以翻转后是 [a0, a_{n-1}, ..., 1]。 % 为了符合公式习惯，我们定义 coeffs(1,:) = [a0, a1, ..., a_{n-1}, 1] (其中a_n=1) % 即 coeffs(m, k) 对应 a_{k-1}^{(m)}， k从1到n+1。 coeffs = zeros(max_iter+1, n+1); coeffs(1, :) = a_prev; % 第1行是初始迭代（m=0） % 3. 执行Graeffe迭代 for m = 1:max_iter a_current = zeros(1, n+1); a_current(end) = 1; % 最高次项系数始终为1 for k = 0:n-1 % k对应公式中的索引，计算 a_k^(m) idx = k + 1; % MATLAB索引 sum_term = 0; for j = 1:min(k, n-k) if (k-j >= 0) && (k+j <= n) prod_term = coeffs(m, k-j+1) * coeffs(m, k+j+1); sum_term = sum_term + ((-1)^j) * prod_term; end end a_current(idx) = ((-1)^k) * (coeffs(m, idx)^2 + 2 * sum_term); end coeffs(m+1, :) = a_current; end这段代码严格实现了递推公式。内层循环的边界条件min(k, n-k)和索引检查(k-j >= 0) && (k+j <= n)确保了不会访问数组越界。这是实现中的第一个关键细节。
3.2 从迭代系数中恢复根的模长
迭代完成后，coeffs矩阵的每一行对应一次迭代的系数。现在，我们需要利用最后两次迭代（s和s-1）的系数来恢复根的模长。
% 4. 从最后两次迭代恢复根的模长 s = max_iter; mod_roots = zeros(1, n); for i = 1:n % 注意公式：|r_i| ≈ ( |a_{n-i}^{(s)} / a_{n-i}^{(s-1)}| )^(1/2^s) % 在我们的coeffs数组中，a_{n-i}^{(s)} 对应 coeffs(s+1, n-i+1) % 因为coeffs行索引从1开始（m=0），所以第s次迭代在 s+1 行。 num = coeffs(s+1, n-i+1); % a_{n-i}^{(s)} den = coeffs(s, n-i+1); % a_{n-i}^{(s-1)} if den ~= 0 ratio = abs(num / den); mod_roots(i) = ratio^(1/(2^s)); else % 如果分母为零，说明对应根的模长可能非常小或计算溢出，需要特殊处理 warning('分母为零在恢复根模长时发生，根索引 %d。这可能意味着该根模长为零或数值下溢。', i); mod_roots(i) = 0; end end % 对模长进行排序（通常从大到小），符合我们恢复的顺序假设 mod_roots = sort(mod_roots, 'descend');这里有几个非常重要的实操心得：
- 分母为零的处理：在迭代后期，某些系数可能变得极其微小甚至下溢为零。这通常意味着对应根的模长非常小（接近零）。直接忽略或赋值为零是一个合理的策略，但最好记录警告，以便后期验证。
- 排序的假设：我们默认
|r_1| > |r_2| > ... > |r_n|，这样系数比a_{n-i}^{(s)} / a_{n-i}^{(s-1)}才近似对应第i大模长的根。如果根模长非常接近，这个对应关系会失效，这是Graeffe方法的主要局限之一。排序操作是为了让输出结果看起来更有序，但并不能解决根本的“根簇”问题。 - 指数运算的数值稳定性：
ratio^(1/(2^s))在s较大（如>20）且ratio极大或极小时，可能引发浮点溢出或下溢。在实际代码中，可以考虑在迭代过程中定期对系数进行缩放（例如，除以某个范数），或者使用对数形式计算：mod = exp( (log(abs(num)) - log(abs(den))) / (2^s) )。
3.3 恢复根的辐角（符号与相位）
仅知道模长还不够，我们需要完整的复数根。对于实系数多项式，根要么是实数，要么是共轭复数对。恢复辐角θ_i（使得r_i = |r_i| * exp(1i*θ_i)）更为微妙。
一种经典方法是利用两次迭代的系数比值来估算辐角的两倍关系，然后通过反推得到。但对于实系数多项式，一个更简单实用的策略是：先假设所有根都是实数，通过符号确定正负，然后用牛顿迭代法或直接代入原多项式进行精修。
% 5. 尝试恢复根的符号（针对实根）并初步构建根 roots_approx = mod_roots; % 先假设都是正实数 % 我们可以利用初始多项式的系数符号关系（笛卡尔符号法则）来辅助判断，但更直接的是用一次代值。 % 一个粗糙但常有效的方法：检查 P(|r|) 和 P(-|r|) 的符号。 for i = 1:n r_pos = roots_approx(i); r_neg = -roots_approx(i); % 计算多项式在正负模长处值 val_pos = polyval(p, r_pos); % p是降幂排列的原始系数 val_neg = polyval(p, r_neg); % 如果改变符号后，多项式的绝对值显著变小，则认为符号取反更优 if abs(val_neg) < abs(val_pos) roots_approx(i) = r_neg; end % 注意：这只对实根有效。对于复根，这样会得到错误结果。 end % 6. 处理共轭复根 % 经过上述步骤，我们得到一组实数候选根。但原多项式可能有复根。 % 一个检查方法是：计算这些候选根构成的多项式，并与原多项式比较残差。 % 如果残差大，说明存在复根。一个更系统的方法是使用恢复的模长和初步的实根猜测， % 然后调用一次MATLAB内置的 `roots` 函数作为精修起点，但这失去了Graeffe的意义。 % 另一种Graeffe特有的方法：利用系数比值恢复辐角。 % 这里展示一个简化的思路：对于疑似复共轭对，其模长应近似相等。 % 我们可以对模长数组进行配对。 tol = 1e-6; % 模长相等容忍度 complex_pairs = []; i = 1; while i < n if abs(mod_roots(i) - mod_roots(i+1)) < tol * max(mod_roots(i), mod_roots(i+1)) % 找到一对模长相等的根，假设是共轭复根 % 需要确定辐角。一个近似公式：θ ≈ (1/2^s) * arg( a_{n-i}^{(s)} / a_{n-i}^{(s-1)} ) % 注意这里索引要对应好 idx_coeff = n - i + 1; num = coeffs(s+1, idx_coeff); den = coeffs(s, idx_coeff); phase = angle(num / den) / (2^s); % 得到的是 2θ？需要仔细推导公式。 % 实际上，对于共轭复根 r 和 r*，其平方后辐角加倍。经过s次迭代后， % 其2^s次幂的辐角是 2^s * θ。而系数比值的辐角包含了这个信息。 % 更稳妥的做法：承认这一步的复杂性，对于复根，我们输出模长和初步相位后， % 强烈建议用牛顿法或内置函数精修。 theta = phase; % 这里是一个示意，具体符号和因子需要根据公式推导调整 r_mod = mod_roots(i); roots_approx(i) = r_mod * exp(1i * theta); roots_approx(i+1) = r_mod * exp(-1i * theta); i = i + 2; % 跳过一对 else i = i + 1; end end我必须坦诚地告诉你，恢复辐角是Graeffe方法实现中最棘手、数值最不稳定的部分。上面的代码只是一个示意框架。在历史文献和实际应用中，对于复根的处理有更复杂的公式（涉及多个系数的比值和符号）。因此，一个更工程化的做法是：
将Graeffe方法视为一个强大的“根定位器”或“初始猜测生成器”。我们利用它稳定地计算出根的模长|r_i|，并大致判断是实根还是复根（通过观察系数比值的虚部或模长配对）。然后，以这些信息作为初始值，调用一个局部收敛的迭代法（如牛顿-拉弗森法、拉盖尔法）进行精修。这样结合了Graeffe的全局稳定性和局部方法的高精度。
3.4 最终的精修与验证步骤
% 7. 使用牛顿法进行精修（可选但推荐） % 使用Graeffe的结果作为牛顿法的初始值 options = optimset('Display', 'off', 'TolFun', 1e-12, 'TolX', 1e-12); for i = 1:n try % 使用fzero处理实根，但需要函数句柄。更通用的是用自定义的牛顿迭代。 % 这里简单演示使用MATLAB的`roots`函数作为“裁判”来对比，实际中应实现牛顿迭代。 % 我们实现一个简单的牛顿迭代： root_guess = roots_approx(i); for iter = 1:20 [P_val, P_der] = polyval_der(p, root_guess); % 需要计算多项式和导数值 if abs(P_der) < 1e-14 break; % 导数太小，避免除零 end step = P_val / P_der; root_guess = root_guess - step; if abs(step) < 1e-12 break; end end roots_approx(i) = root_guess; catch % 如果精修失败，保留Graeffe的结果 end end % 8. 排序并返回结果（通常按实部或模长排序） [~, idx] = sort(real(roots_approx)); roots_approx = roots_approx(idx); end % 辅助函数：计算多项式在某点的值和导数值 function [val, der] = polyval_der(coeff, x) % coeff: 降幂排列，如 [1, -5, 6] n = length(coeff) - 1; val = coeff(1); der = 0; for i = 1:n der = der * x + val; val = val * x + coeff(i+1); end end4. 实战测试、对比分析与经典“踩坑”点
现在，让我们用几个例子来测试我们的graeffeRoots函数，并与MATLAB内置的roots函数对比，同时指出过程中会遇到哪些坑。
4.1 测试案例一：良态多项式
% 案例1：根为 2, 3 的简单二次方程 p1 = [1, -5, 6]; % x^2 -5x +6 r_graeffe1 = graeffeRoots(p1, 8); r_matlab1 = roots(p1); disp('案例1: x^2 -5x +6 = 0'); disp(['Graeffe结果: ', num2str(r_graeffe1.')]); disp(['MATLAB roots结果: ', num2str(r_matlab1.')]); disp(['误差范数: ', num2str(norm(sort(r_graeffe1)-sort(r_matlab1)))]);对于这种简单情况，只要迭代次数足够（比如8次），Graeffe方法应该能给出相当精确的结果，误差在1e-10量级或更低。精修步骤会进一步将误差降到机器精度附近。
4.2 测试案例二：具有大动态范围根的多项式
这是Graeffe方法最能发挥优势的场景。
% 案例2：根模长差异巨大，例如 (x-10)(x-0.1)(x+0.01) = x^3 - 10.09x^2 + 1.011x - 0.001 p2 = [1, -10.09, 1.011, -0.001]; r_graeffe2 = graeffeRoots(p2, 10); r_matlab2 = roots(p2); disp('案例2: 根模长差异大 (10, 0.1, 0.01)'); disp(['Graeffe结果: ', num2str(r_graeffe2.')]); disp(['MATLAB roots结果: ', num2str(r_matlab2.')]);你会发现，对于这个多项式，roots函数可能在小根（0.01）上出现相对较大的误差，因为大根（10）在计算中可能引起数值抵消。而Graeffe方法通过根平方法放大了差距，反而能更稳定地捕捉到小根。这是一个关键洞察：对于病态多项式（如根尺度差异大），基于伴侣矩阵特征值的roots函数可能不如Graeffe方法稳健。
4.3 测试案例三：重根或接近的根（Graeffe的“天敌”）
% 案例3：接近的根或重根，例如 (x-1)^2 * (x-2) = x^3 -4x^2 +5x -2 p3 = [1, -4, 5, -2]; r_graeffe3 = graeffeRoots(p3, 8); r_matlab3 = roots(p3); disp('案例3: 重根 (1, 1, 2)'); disp(['Graeffe结果: ', num2str(r_graeffe3.')]); disp(['MATLAB roots结果: ', num2str(r_matlab3.')]);这里就是大坑所在！Graeffe方法的基本假设是根的模长在经过多次平方后能被充分分离。如果两个根模长相等（如重根|r1|=|r2|）或者非常接近，这个假设就不成立了。迭代后，对应的系数比值将不收敛到单个值，而是出现振荡或无法给出清晰信息。我们的恢复公式将失效，计算出的模长可能完全错误。对于重根(x-1)^2，两个根都是1，模长相等，Graeffe方法会遭遇严重困难。此时，算法可能输出两个不相等且偏离1的模长值。
踩坑心得与应对策略：
- 检测根簇：在迭代过程中，可以监控相邻恢复的模长比值。如果连续几次迭代，某两个模长的比值始终接近1，那么它们很可能属于一个根簇（重根或非常接近的根）。
- 提前终止或切换策略：一旦检测到根簇，应提前终止Graeffe迭代。对于根簇，更好的方法是使用专门处理重根的算法，或者将多项式除以已求得的单根因子后降阶，再对降阶后的多项式（其根簇可能被分离）应用Graeffe或其他方法。
- 不要迷信单一方法：Graeffe是一个强大的工具，但绝非万能。在实际应用中，应将它与
roots、fzero（针对实根）等函数结合使用。例如，可以先尝试roots，如果结果数值表现很差（如残差很大），再换用Graeffe求取初始猜测，然后用牛顿法精修。
4.4 测试案例四：高次多项式与数值溢出
% 案例4：高次多项式，系数范围大 r_large = [100, -50, 3, 0.01, -0.005]; p4 = poly(r_large); % poly函数根据根生成多项式系数 r_graeffe4 = graeffeRoots(p4, 12); r_matlab4 = roots(p4); disp('案例4: 高次多项式，根动态范围大'); % 比较残差：计算多项式在求得根处的值 residual_graeffe = max(abs(polyval(p4, r_graeffe4))); residual_matlab = max(abs(polyval(p4, r_matlab4))); disp(['Graeffe最大残差: ', num2str(residual_graeffe)]); disp(['MATLAB roots最大残差: ', num2str(residual_matlab)]);对于高次多项式，Graeffe迭代过程中的系数可能会发生数量级的剧烈变化，导致浮点数上溢（变成Inf）或下溢（变成0）。这是我们代码中另一个潜在的坑。
避坑技巧：系数缩放在每次迭代后，可以对系数向量进行缩放，例如除以该次迭代系数的最大绝对值（或某种范数）。这不会改变根的比例关系（因为整个多项式被缩放了一个常数因子），但能有效将系数控制在合理的范围内，避免数值溢出。
% 在Graeffe迭代循环内，计算完a_current后添加： scale_factor = max(abs(a_current)); if scale_factor > 1e100 || scale_factor < 1e-100 % 进行缩放 a_current = a_current / scale_factor; % 注意：缩放会影响后续系数比值的计算，需要在恢复根模长时补偿这个缩放因子。 % 一种方法是记录每次迭代的缩放因子 log_scale(m) = log(scale_factor)。 % 恢复模长时，公式变为 |r_i| ≈ exp( (log|num| - log|den| + log_scale(s) - log_scale(s-1)) / 2^s ) end实现完整的缩放逻辑需要小心处理，但它能极大增强算法对高次多项式的鲁棒性。
5. 总结与工程应用建议
经过上面的原理剖析、代码实现和测试分析，你应该对Graeffe根平方法有了一个从理论到实践的全面认识。它不是一个“黑盒”函数，而是一种蕴含着经典数值智慧的思想。在现代计算环境中，它的角色更像是一个“特种兵”，在特定战场（病态多项式、大动态范围根）上表现卓越，而非日常的“通用士兵”。
给实际工程应用的建议：
定位，而非终结：将Graeffe方法视为获取高质量初始猜测的工具。用它的输出（尤其是模长信息）作为牛顿法、拉盖尔法或
roots函数（在降阶后）的起点，往往能解决直接用这些方法收敛失败或不稳定的问题。迭代次数的选择：迭代次数
s并非越多越好。通常8-12次足以将根的模长差异放大到足够大。过多迭代会增加数值溢出风险和计算量。一个实用的策略是监控系数比值的变化，当相邻迭代恢复的模长变化小于某个阈值时停止。复根处理的务实态度：如前所述，完整、稳定地恢复复根辐角是复杂的。对于必须处理大量复根的场景，依赖Graeffe完成全部工作可能不划算。更常见的做法是：先用Graeffe判断是否存在明显的主导实根或模长差异大的根，将它们提取出来（通过多项式降阶），对剩下的多项式（其次数降低，病态程度可能减轻）再使用
roots函数。与MATLAB生态结合：在MATLAB中，你可以将
graeffeRoots函数封装成一个接收p和可选参数（最大迭代次数、容忍度等）的函数。在函数内部实现我们讨论的缩放、根簇检测、以及可选的牛顿精修。这样，你就拥有了一个属于自己的、增强版的求根工具库。
最后，我想分享一点个人体会：学习像Graeffe这样的经典算法，最大的收获往往不是算法本身，而是那种“将困难问题转化为可解问题”的思维模式。在面对一个数值不稳定的系统时，我们是否可以像Graeffe那样，通过某种变换改变问题的“尺度”或“视角”，从而暴露出隐藏的结构？这种思维，在调试一个发散的优化算法、分析一个震荡的控制系统时，同样珍贵。希望这篇文章不仅给了你一段可运行的代码，更提供了一种解决数值问题的思路。