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^(m1) (-1)^k * [ (a_k^(m))^2 2 * Σ_{j1}^{min(k, n-k)} (-1)^j * a_{k-j}^(m) * a_{kj}^(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_n1) % 即 coeffs(m, k) 对应 a_{k-1}^{(m)} k从1到n1。 coeffs zeros(max_iter1, n1); coeffs(1, :) a_prev; % 第1行是初始迭代m0 % 3. 执行Graeffe迭代 for m 1:max_iter a_current zeros(1, n1); 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) (kj n) prod_term coeffs(m, k-j1) * coeffs(m, kj1); sum_term sum_term ((-1)^j) * prod_term; end end a_current(idx) ((-1)^k) * (coeffs(m, idx)^2 2 * sum_term); end coeffs(m1, :) a_current; end这段代码严格实现了递推公式。内层循环的边界条件min(k, n-k)和索引检查(k-j 0) (kj 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(s1, n-i1) % 因为coeffs行索引从1开始m0所以第s次迭代在 s1 行。 num coeffs(s1, n-i1); % a_{n-i}^{(s)} den coeffs(s, n-i1); % 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(i1)) tol * max(mod_roots(i), mod_roots(i1)) % 找到一对模长相等的根假设是共轭复根 % 需要确定辐角。一个近似公式θ ≈ (1/2^s) * arg( a_{n-i}^{(s)} / a_{n-i}^{(s-1)} ) % 注意这里索引要对应好 idx_coeff n - i 1; num coeffs(s1, 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(i1) 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(i1); 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)(x0.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那样通过某种变换改变问题的“尺度”或“视角”从而暴露出隐藏的结构这种思维在调试一个发散的优化算法、分析一个震荡的控制系统时同样珍贵。希望这篇文章不仅给了你一段可运行的代码更提供了一种解决数值问题的思路。