61次多项式在sqrt(2)附近求值出现数值异常的原因及解决方法

问题描述 #

使用Python计算由P.T.P Tang推广的Remez算法生成的61次非经典Chebyshev多项式，该多项式目标是在复平面紧子集上最小化||zⁿ+低阶项||。但在√2附近的虚数方向（具体为I = np.sqrt(2)+1j*np.linspace(0,1e-6,1000, dtype = np.clongdouble)）求值时，多项式绝对值的图像出现异常波动。

复现代码 #

import numpy as np
import matplotlib.pyplot as plt

I = np.sqrt(2)+1j*np.linspace(0,1e-6,1000, dtype = np.clongdouble)
Y = np.abs(np.polynomial.polynomial.polyval(I, coefficients_rev))
plt.plot(np.linspace(0,1,1000),Y)
plt.show()

多项式系数 #

coefficients = np.array([ 1.00000000e+00+0.j,  0.00000000e+00+0.j, -3.04750157e+01+0.j,
    0.00000000e+00+0.j,  4.49106800e+02+0.j,  0.00000000e+00+0.j,
   -4.26240239e+03+0.j,  0.00000000e+00+0.j,  2.92734180e+04+0.j,
    0.00000000e+00+0.j, -1.54973591e+05+0.j,  0.00000000e+00+0.j,
    6.57829830e+05+0.j,  0.00000000e+00+0.j, -2.29933795e+06+0.j,
    0.00000000e+00+0.j,  6.74451933e+06+0.j,  0.00000000e+00+0.j,
   -1.68346568e+07+0.j,  0.00000000e+00+0.j,  3.61319542e+07+0.j,
    0.00000000e+00+0.j, -6.72090548e+07+0.j,  0.00000000e+00+0.j,
    1.08986264e+08+0.j,  0.00000000e+00+0.j, -1.54736506e+08+0.j,
    0.00000000e+00+0.j,  1.92921639e+08+0.j,  0.00000000e+00+0.j,
   -2.11600531e+08+0.j,  0.00000000e+00+0.j,  2.04319983e+08+0.j,
    0.00000000e+00+0.j, -1.73627642e+08+0.j,  0.00000000e+00+0.j,
    1.29668219e+08+0.j,  0.00000000e+00+0.j, -8.48892421e+07+0.j,
    0.00000000e+00+0.j,  4.85309319e+07+0.j,  0.00000000e+00+0.j,
   -2.41002152e+07+0.j,  0.00000000e+00+0.j,  1.03215590e+07+0.j,
    0.00000000e+00+0.j, -3.77608656e+06+0.j,  0.00000000e+00+0.j,
    1.16509284e+06+0.j,  0.00000000e+00+0.j, -2.97955777e+05+0.j,
    0.00000000e+00+0.j,  6.16343028e+04+0.j,  0.00000000e+00+0.j,
   -9.94806340e+03+0.j,  0.00000000e+00+0.j,  1.18267487e+03+0.j,
    0.00000000e+00+0.j, -9.31129650e+01+0.j,  0.00000000e+00+0.j,
    3.72597504e+00+0.j,  0.00000000e+00+0.j], dtype=np.clongdouble)

coefficients_rev = coefficients[::-1]

数值异常成因分析 #

• 高阶单项式基的数值不稳定性：61次多项式的系数量级从1e0到1e8跨度极大，使用单项式基求值时，高阶项的数值会掩盖低阶项；同时√2附近可能是多项式的近根或极值边界，微小的虚部扰动会引发各阶项的相位干涉，放大舍入误差。
• 灾难性抵消与舍入误差：在√2附近，多项式各阶项可能存在相互抵消的情况，即使是clongdouble精度的浮点数也无法准确捕捉这种抵消，导致计算结果出现无规则波动。
• 复平面局部敏感性：该区域属于多项式在复平面上的“陡峭”变化区，函数值对自变量的微小变化极度敏感，有限精度的数值计算无法复现理论上的平滑变化，反而将舍入误差放大为可见波动。

可行解决办法 #

• 转换为切比雪夫基求值：切比雪夫多项式基的数值稳定性远优于单项式基，将原多项式转换为切比雪夫基表示后再求值，能大幅降低系数量级差异带来的误差。
• 实现带误差补偿的霍纳法：手动实现融入Kahan求和的霍纳算法，减少求值过程中的舍入误差累积，比np.polyval的默认实现更稳定。
• 局部低阶近似：在√2的小邻域内对多项式做泰勒展开，用低阶多项式近似原高阶多项式，减少求值项数和误差来源。
• 因式分解优化：若√2是多项式的近根，可将多项式分解为(z - √2) * Q(z)的形式，先计算Q(z)的绝对值再乘以|z - √2|，避免直接计算高阶项的抵消。
• 使用任意精度计算库：借助mpmath等支持任意精度的库，提高计算精度，压制舍入误差的影响。

内容的提问来源于stack exchange，提问作者Olof R