NumPy Polynomial 便捷类完全指南从幂级数到 Chebyshev 拟合的实战手册【免费下载链接】numpyThe fundamental package for scientific computing with Python.项目地址: https://gitcode.com/gh_mirrors/nu/numpynumpy.polynomial子包仓库中位于 numpy/polynomial为科学计算提供了六种开箱即用的多项式便捷类convenience classes它们共享同一套面向对象的 APIPolynomial、Chebyshev、Legendre、Laguerre、Hermite与HermiteE。本指南以官方参考文档 routines.polynomials.classes.rst 为骨架结合仓库源码深入讲解这些类的构造、算术运算、微积分、类型转换与最小二乘拟合并重点剖析domain/window映射机制这一拟合成败的关键。读完本文你将能够熟练使用统一 API 完成多项式建模、求根、拟合与序列转换。六种便捷类一览统一的序列模型多项式便捷类分别位于numpy.polynomial包下的六个子模块中每个类代表一种正交多项式基函数展开的有限级数类名提供内容源码位置Polynomial幂级数Power seriesnumpy/polynomial/polynomial.pyChebyshevChebyshev 级数numpy/polynomial/chebyshev.pyLegendreLegendre 级数numpy/polynomial/legendre.pyLaguerreLaguerre 级数numpy/polynomial/laguerre.pyHermiteHermite 级数物理学家版本numpy/polynomial/hermite.pyHermiteEHermiteE 级数概率论版本numpy/polynomial/hermite_e.py这里的“级数”指对应基函数乘以系数后的有限和。例如幂级数$$ p(x) 1 2x 3x^2 $$其系数为[1, 2, 3]。同样的系数写成 Chebyshev 级数则是$$ p(x) 1,T_0(x) 2,T_1(x) 3,T_2(x) $$更一般地$$ p(x) \sum_{i0}^{n} c_i,T_i(x) $$其中 $T_n$ 是 $n$ 次 Chebyshev 函数但也可以是其他任何类的基函数。所有类的统一约定是系数c[i]与 $i$ 次基函数对应参见 chebyshev.py 中基类命名basis_name T。所有类都不可变immutable且拥有相同的方法集合。特别是它们实现了 Python 数值运算符、-、*、//、%、divmod、**、和!——后两者因浮点舍入误差可能有些棘手。统一的抽象基类 ABCPolyBase六种类的共同 API 全部由抽象基类 ABCPolyBase 提供。它声明了每个具体类必须实现的静态方法_add、_sub、_mul、_div、_pow、_val、_int、_der、_fit、_line、_roots、_fromroots以及类属性domain、window、basis_name、maxpower。从源码结构可以推断各子模块只需实现这些与基函数相关的底层运算就能自动获得全部高层功能——这正是六类 API 完全一致的架构原因。基础操作构造、打印与属性构造实例与三大属性可以从numpy.polynomial包直接导入类也可以从对应类型的子模块导入。官方文档以最熟悉的Polynomial为例 from numpy.polynomial import Polynomial as P p P([1,2,3]) p Polynomial([1., 2., 3.], domain[-1., 1.], window[-1., 1.], symbolx)长格式打印包含三部分系数coefficients、域domain、窗window新版还包含符号symbol p.coef array([1., 2., 3.]) p.domain array([-1., 1.]) p.window array([-1., 1.])各属性的语义在 ABCPolyBase 的文档字符串 中有精确定义coef按次数递增排列的系数数组(1, 2, 3)表示1*P_0(x) 2*P_1(x) 3*P_2(x)domain数据所在区间[domain[0], domain[1]]通过平移缩放线性映射到windowwindow映射目标区间即多项式表现良好的区间symbol字符串表示中自变量符号必须是合法的 Python 标识符默认x1.24 版本新增见 构造器校验逻辑。字符串表示与打印风格切换print(p)输出更熟悉的数学表达式 print(p) 1.0 2.0·x 3.0·x²字符串表示默认使用 Unicode 字符表达幂次与下标Windows 除外因部分 Windows 终端默认字体不支持完整 Unicode 上下标范围见 _polybase.py 中_use_unicode的平台判断。ASCII 表示法同样可用并可通过包级函数set_default_printstyle全局切换 np.polynomial.set_default_printstyle(ascii) print(p) 1.0 2.0 x 3.0 x**2该函数实现在 numpy/polynomial/init.py仅接受ascii与unicode两个取值传入其他值会抛出ValueError。也可对单个实例通过字符串格式化控制且格式化优先级高于全局默认 print(f{p:unicode}) 1.0 2.0·x 3.0·x²代数与算术运算domain与window暂时忽略先看各类基础代数运算以下结果均来自官方文档演示。加法与减法 p p Polynomial([2., 4., 6.], domain[-1., 1.], window[-1., 1.], symbolx) p - p Polynomial([0.], domain[-1., 1.], window[-1., 1.], symbolx)注意p - p的系数被压缩为[0.]说明减法结果自动修剪了尾随零。乘法与幂 p * p Polynomial([ 1., 4., 10., 12., 9.], domain[-1., 1.], window[-1., 1.], symbolx) p**2 Polynomial([ 1., 4., 10., 12., 9.], domain[-1., 1.], window[-1., 1.], symbolx)从基类实现看这些运算都经由_get_coefficients提取运算对象系数后调用对应的_add/_mul/_pow完成见 _polybase.py 的__add__/__mul__/__pow__。此外基类还定义了maxpower 100来限制幂运算导致的级数规模失控T_n^m的次数为n*m见 _polybase.py#L76-L77。除法、取余与 divmod整除//是多项式类的除法运算符多项式在此被当作类似整数的对象处理。Python 3.x 中/仅对除以标量有效未来会弃用而//、%、divmod均可用 p // P([-1, 1]) Polynomial([5., 3.], domain[-1., 1.], window[-1., 1.], symbolx) p % P([-1, 1]) Polynomial([6.], domain[-1., 1.], window[-1., 1.], symbolx) quo, rem divmod(p, P([-1, 1])) quo Polynomial([5., 3.], domain[-1., 1.], window[-1., 1.], symbolx) rem Polynomial([6.], domain[-1., 1.], window[-1., 1.], symbolx)__floordiv__与__mod__在源码中都是委托给__divmod__完成的见 _polybase.py#L565-L575__divmod__内部调用self._div得到商与余数见 _polybase.py#L577-L587。求值与替换函数复合多项式实例可以直接当函数调用对数组逐元素求值 x np.arange(5) p(x) array([ 1., 6., 17., 34., 57.]) x np.arange(6).reshape(3,2) p(x) array([[ 1., 6.], [17., 34.], [57., 86.]])替换substitution把多项式代入x并展开若把多项式视为函数这就是函数复合 p(p) Polynomial([ 6., 16., 36., 36., 27.], domain[-1., 1.], window[-1., 1.], symbolx)求根 p.roots() array([-0.33333333-0.47140452j, -0.333333330.47140452j])roots()返回复根数组。从源码看它调用各类型的_roots后通过pu.mapdomain把根从window映射回domain见 _polybase.py#L900-L913文档字符串也提醒根越偏离domain区间精度越低。隐式类型转换与混合运算规则不必总是显式构造Polynomial实例元组、列表、数组和标量会在算术运算中被自动转换为系数参与计算 p [1, 2, 3] Polynomial([2., 4., 6.], domain[-1., 1.], window[-1., 1.], symbolx) [1, 2, 3] * p Polynomial([ 1., 4., 10., 12., 9.], domain[-1., 1.], window[-1., 1.], symbolx) p / 2 Polynomial([0.5, 1. , 1.5], domain[-1., 1.], window[-1., 1.], symbolx)domain、window 或类型不同的多项式不能直接混合运算 from numpy.polynomial import Chebyshev as T p P([1], domain[0,1]) TypeError: Domains differ p P([1], window[0,1]) TypeError: Windows differ p T([1]) TypeError: Polynomial types differ这些异常正是由_get_coefficients中的兼容性检查抛出的见 _polybase.py#L256-L290先检查类型是否一致再检查domain、window、symbol是否分别相同。利用替换实现类间转换不同类型不能直接做算术但可以用作替换substitution实现转换。事实上convert正是通过“把kind类的恒等多项式代入自身”来实现类型、domain、window 三重转换的见 convert 实现 p(T([0, 1])) Chebyshev([2.5, 2. , 1.5], domain[-1., 1.], window[-1., 1.], symbolx)这得到多项式p的 Chebyshev 形式。其原理是 $T_1(x) x$代入x不改变原多项式但所有乘除都在 Chebyshev 级数下进行因此结果类型为Chebyshev。不可变性与增强赋值所有多项式实例设计为不可变因此、-等增强赋值操作以及其他会破坏不可变性的功能都被有意未实现unimplemented。基类文档也声明__hash__ None并关闭了 ndarray 子类的 ufunc 与运算符交互见 _polybase.py#L70-L74。微积分微分与积分多项式实例支持积分与微分注意积分默认的下限与积分常数均为 0但都可自定义 from numpy.polynomial import Polynomial as P p P([2, 6]) p.integ() # 积分一次 Polynomial([0., 2., 3.], domain[-1., 1.], window[-1., 1.], symbolx) p.integ(2) # 积分两次 Polynomial([0., 0., 1., 1.], domain[-1., 1.], window[-1., 1.], symbolx) p.integ(lbnd-1) # 指定积分下限 Polynomial([-1., 2., 3.], domain[-1., 1.], window[-1., 1.], symbolx) p.integ(lbnd-1, k1) # 下限为 -1 且积分常数为 1 Polynomial([0., 2., 3.], domain[-1., 1.], window[-1., 1.], symbolx)m积分次数非负整数k积分常数数组k[0]应用于第一次积分、k[1]应用于第二次依此类推缺失值补零见 integ 签名。微分更简单唯一参数是微分次数 p P([1, 2, 3]) p.deriv(1) Polynomial([2., 6.], domain[-1., 1.], window[-1., 1.], symbolx) p.deriv(2) Polynomial([6.], domain[-1., 1.], window[-1., 1.], symbolx)值得一提的实现细节integ与deriv都先通过mapparms()获取domain→window的线性映射参数(off, scl)积分下限与微分结果都会按映射关系做缩放1./scl保证结果仍在该实例的domain上语义正确见 _polybase.py#L870-L898。其他构造方式fromroots、convert、basis 与 cast按系数构造只是获得多项式实例的一种方式还可以按根构造、从其他类型转换、以及通过最小二乘拟合构造拟合单独成节。按根构造与类型转换 from numpy.polynomial import Polynomial as P from numpy.polynomial import Chebyshev as T p P.fromroots([1, 2, 3]) p Polynomial([-6., 11., -6., 1.], domain[-1., 1.], window[-1., 1.], symbolx) p.convert(kindT) Chebyshev([-9. , 11.75, -3. , 0.25], domain[-1., 1.], window[-1., 1.], symbolx)fromroots构造的是 $\prod_i (x - r_i)$ 的展开其内部会先把根从domain映射到window再生成系数见 _polybase.py#L1037-L1077。convert方法还能同时转换 domain 与 window p.convert(kindT, domain[0, 1]) Chebyshev([-2.4375 , 2.96875, -0.5625 , 0.03125], domain[0., 1.], window[-1., 1.], symbolx) p.convert(kindP, domain[0, 1]) Polynomial([-1.875, 2.875, -1.125, 0.125], domain[0., 1.], window[-1., 1.], symbolx)basis 与 cast 类方法NumPy 1.7.0 起新增basis与cast类方法。cast与convert类似basis返回指定次数的基多项式 P.basis(3) Polynomial([0., 0., 0., 1.], domain[-1., 1.], window[-1., 1.], symbolx) T.cast(p) Chebyshev([-9. , 11.75, -3. , 0.25], domain[-1., 1.], window[-1., 1.], symbolx)basis的实现是把deg位置置 1、其余置 0 的系数交给构造器见 _polybase.py#L1114-L1151cast则委托给实例方法convert见 _polybase.py#L1153-L1191。类型转换的精度警告类型之间的转换虽然有用但不建议常规使用把 50 次的 Chebyshev 级数转换为同次幂级数时数值精度损失可能使求值结果近乎随机。基类的convert文档也明确提示“域和类型之间的转换可能产生数值上病态的级数”见 _polybase.py#L802-L806。拟合Fittingdomain 与 window 的用武之地拟合正是domain与window两个属性存在的根本原因。问题背景如下$[-1, 1]$ 区间内前几阶 Chebyshev 多项式是位于 ±1 之间的等波纹equiripple函数形态良好但同样的曲线放到 $[-2, 2]$ 区间良好区域就缩小到几乎可忽略。机制用 Chebyshev 多项式拟合时我们希望利用 $x$ 位于 $[-1, 1]$ 的“良好区域”——这正是window指定的区间。但待拟合数据几乎不可能恰好全部落在这个区间内于是用domain指定数据点所在的区间。拟合时先将domain通过线性变换映射到window再对映射后的数据点做常规最小二乘拟合。拟合结果的window与domain会随返回的级数一起保存并在后续求值、微分等操作中自动生效。若调用时未指定拟合例程会使用默认window和覆盖全部数据点的最小domain。线性映射的具体实现映射的数学核心在 polyutils.py 的三个函数中getdomain(x)返回覆盖数据的最小区间实数据为[min, max]复数据为包含全部点、边与坐标轴对齐的最小矩形的左下/右上角见 polyutils.py#L194-L239mapparms(old, new)返回线性映射L(x) off scl*x的参数满足old[i] → new[i]见 polyutils.py#L241-L286mapdomain(x, old, new)把映射应用到数据点x见 polyutils.py#L288-L323。fit的流程在 _polybase.py#L1014-L1034 中清晰可见domain缺省时由pu.getdomain(x)自动推导若数据退化为单点则自动扩展 ±1然后xnew pu.mapdomain(x, domain, window)映射数据点再交给具体类的_fit完成最小二乘。对含噪正弦数据的拟合示例 import numpy as np import matplotlib.pyplot as plt from numpy.polynomial import Chebyshev as T np.random.seed(11) x np.linspace(0, 2*np.pi, 20) y np.sin(x) np.random.normal(scale.1, sizex.shape) p T.fit(x, y, 5) plt.plot(x, y, o) xx, yy p.linspace() plt.plot(xx, yy, lw2) p.domain array([0. , 6.28318531]) p.window array([-1., 1.]) plt.show()注意因为调用时未显式指定domain拟合出的p.domain自动取为覆盖数据的最小区间[0, 2π]而p.window保持类默认的[-1, 1]——两者正是通过getdomain与mapparms联系起来的。fit 的完整参数fit类方法的签名见 _polybase.py#L946-L1012比文档示例更丰富x, y形如(M,)的样本点坐标deg拟合多项式的次数int表示包含直到deg次的所有项NumPy ≥ 1.11 还支持传入整数列表来指定要包含的特定次数项domainNone时自动选覆盖x的最小domain[]时用类默认domain[]选项 1.5.0 新增rcond相对条件数小于该值相对最大奇异值的奇异值被忽略默认len(x)*eps约 2e-16 量级fullTrue时额外返回[resid, rank, sv, rcond]诊断信息残差平方和、缩放 Vandermonde 矩阵的数值秩、奇异值、rcond详见linalg.lstsqw权重数组w[i]作用于x[i]处的未平方残差y[i] - y_hat[i]逆方差加权时取w[i] 1/sigma(y[i])window返回级数的window默认类默认windowsymbol自变量符号默认x。若需要未缩放、未平移的基多项式系数可执行new_series.convert().coef。绘图的得力助手linspace拟合示例中的p.linspace()是linspace(n100, domainNone)方法返回domain默认用实例自身的domain上n个等距点的(x, y)对主要用作绘图辅助见 _polybase.py#L915-L943。实用补充方法除本文已展开的方法外ABCPolyBase还提供一批常用工具方法见基类degree()返回次数系数长度减 1不检查非零尾项配合trim()可得到真实次数trim(tol0)移除绝对值不超过tol的尾随系数返回新实例cutdeg(deg)/truncate(size)按次数/长度截断级数常用于最小二乘中丢弃高阶小系数copy()返回副本mapparms()返回当前domain→window的线性映射参数(off, scl)has_samecoef/has_samedomain/has_samewindow/has_sametype一致性检查对应算术运算前的兼容性判定identity()恒等多项式p(x) x也是convert的实现基础。这些能力都有对应的测试覆盖例如 numpy/polynomial/tests/test_classes.py 与 test_chebyshev.py 等按类型划分的测试文件可作为进一步研读各类型底层_fit、_roots等实现行为的入口。小结NumPy 多项式便捷类以ABCPolyBase抽象基类统一了六种正交多项式级数的对象模型系数按次数递增存放、实例不可变、完整支持 - * // % divmod **运算与求值/复合/求根并提供integ/deriv微积分、fromroots/convert/basis/cast多种构造路径。真正发挥这些类威力的是domain/window映射机制——它把任意区间上的数据线性映射到基函数表现良好的 $[-1,1]$ 区间再做最小二乘拟合并在返回的实例中持久保存映射关系使后续求值、微分自动一致。理解这一机制是使用Chebyshev.fit等高阶拟合能力时避免数值病态的关键。【免费下载链接】numpyThe fundamental package for scientific computing with Python.项目地址: https://gitcode.com/gh_mirrors/nu/numpy创作声明：本文部分内容由AI辅助生成（AIGC），仅供参考 NumPy Polynomial 便捷类完全指南：从幂级数到 Chebyshev 拟合的实战手册
NumPy Polynomial 便捷类完全指南从幂级数到 Chebyshev 拟合的实战手册【免费下载链接】numpyThe fundamental package for scientific computing with Python.项目地址: https://gitcode.com/gh_mirrors/nu/numpynumpy.polynomial子包仓库中位于 numpy/polynomial为科学计算提供了六种开箱即用的多项式便捷类convenience classes它们共享同一套面向对象的 APIPolynomial、Chebyshev、Legendre、Laguerre、Hermite与HermiteE。本指南以官方参考文档 routines.polynomials.classes.rst 为骨架结合仓库源码深入讲解这些类的构造、算术运算、微积分、类型转换与最小二乘拟合并重点剖析domain/window映射机制这一拟合成败的关键。读完本文你将能够熟练使用统一 API 完成多项式建模、求根、拟合与序列转换。六种便捷类一览统一的序列模型多项式便捷类分别位于numpy.polynomial包下的六个子模块中每个类代表一种正交多项式基函数展开的有限级数类名提供内容源码位置Polynomial幂级数Power seriesnumpy/polynomial/polynomial.pyChebyshevChebyshev 级数numpy/polynomial/chebyshev.pyLegendreLegendre 级数numpy/polynomial/legendre.pyLaguerreLaguerre 级数numpy/polynomial/laguerre.pyHermiteHermite 级数物理学家版本numpy/polynomial/hermite.pyHermiteEHermiteE 级数概率论版本numpy/polynomial/hermite_e.py这里的“级数”指对应基函数乘以系数后的有限和。例如幂级数$$ p(x) 1 2x 3x^2 $$其系数为[1, 2, 3]。同样的系数写成 Chebyshev 级数则是$$ p(x) 1,T_0(x) 2,T_1(x) 3,T_2(x) $$更一般地$$ p(x) \sum_{i0}^{n} c_i,T_i(x) $$其中 $T_n$ 是 $n$ 次 Chebyshev 函数但也可以是其他任何类的基函数。所有类的统一约定是系数c[i]与 $i$ 次基函数对应参见 chebyshev.py 中基类命名basis_name T。所有类都不可变immutable且拥有相同的方法集合。特别是它们实现了 Python 数值运算符、-、*、//、%、divmod、**、和!——后两者因浮点舍入误差可能有些棘手。统一的抽象基类 ABCPolyBase六种类的共同 API 全部由抽象基类 ABCPolyBase 提供。它声明了每个具体类必须实现的静态方法_add、_sub、_mul、_div、_pow、_val、_int、_der、_fit、_line、_roots、_fromroots以及类属性domain、window、basis_name、maxpower。从源码结构可以推断各子模块只需实现这些与基函数相关的底层运算就能自动获得全部高层功能——这正是六类 API 完全一致的架构原因。基础操作构造、打印与属性构造实例与三大属性可以从numpy.polynomial包直接导入类也可以从对应类型的子模块导入。官方文档以最熟悉的Polynomial为例 from numpy.polynomial import Polynomial as P p P([1,2,3]) p Polynomial([1., 2., 3.], domain[-1., 1.], window[-1., 1.], symbolx)长格式打印包含三部分系数coefficients、域domain、窗window新版还包含符号symbol p.coef array([1., 2., 3.]) p.domain array([-1., 1.]) p.window array([-1., 1.])各属性的语义在 ABCPolyBase 的文档字符串 中有精确定义coef按次数递增排列的系数数组(1, 2, 3)表示1*P_0(x) 2*P_1(x) 3*P_2(x)domain数据所在区间[domain[0], domain[1]]通过平移缩放线性映射到windowwindow映射目标区间即多项式表现良好的区间symbol字符串表示中自变量符号必须是合法的 Python 标识符默认x1.24 版本新增见 构造器校验逻辑。字符串表示与打印风格切换print(p)输出更熟悉的数学表达式 print(p) 1.0 2.0·x 3.0·x²字符串表示默认使用 Unicode 字符表达幂次与下标Windows 除外因部分 Windows 终端默认字体不支持完整 Unicode 上下标范围见 _polybase.py 中_use_unicode的平台判断。ASCII 表示法同样可用并可通过包级函数set_default_printstyle全局切换 np.polynomial.set_default_printstyle(ascii) print(p) 1.0 2.0 x 3.0 x**2该函数实现在 numpy/polynomial/init.py仅接受ascii与unicode两个取值传入其他值会抛出ValueError。也可对单个实例通过字符串格式化控制且格式化优先级高于全局默认 print(f{p:unicode}) 1.0 2.0·x 3.0·x²代数与算术运算domain与window暂时忽略先看各类基础代数运算以下结果均来自官方文档演示。加法与减法 p p Polynomial([2., 4., 6.], domain[-1., 1.], window[-1., 1.], symbolx) p - p Polynomial([0.], domain[-1., 1.], window[-1., 1.], symbolx)注意p - p的系数被压缩为[0.]说明减法结果自动修剪了尾随零。乘法与幂 p * p Polynomial([ 1., 4., 10., 12., 9.], domain[-1., 1.], window[-1., 1.], symbolx) p**2 Polynomial([ 1., 4., 10., 12., 9.], domain[-1., 1.], window[-1., 1.], symbolx)从基类实现看这些运算都经由_get_coefficients提取运算对象系数后调用对应的_add/_mul/_pow完成见 _polybase.py 的__add__/__mul__/__pow__。此外基类还定义了maxpower 100来限制幂运算导致的级数规模失控T_n^m的次数为n*m见 _polybase.py#L76-L77。除法、取余与 divmod整除//是多项式类的除法运算符多项式在此被当作类似整数的对象处理。Python 3.x 中/仅对除以标量有效未来会弃用而//、%、divmod均可用 p // P([-1, 1]) Polynomial([5., 3.], domain[-1., 1.], window[-1., 1.], symbolx) p % P([-1, 1]) Polynomial([6.], domain[-1., 1.], window[-1., 1.], symbolx) quo, rem divmod(p, P([-1, 1])) quo Polynomial([5., 3.], domain[-1., 1.], window[-1., 1.], symbolx) rem Polynomial([6.], domain[-1., 1.], window[-1., 1.], symbolx)__floordiv__与__mod__在源码中都是委托给__divmod__完成的见 _polybase.py#L565-L575__divmod__内部调用self._div得到商与余数见 _polybase.py#L577-L587。求值与替换函数复合多项式实例可以直接当函数调用对数组逐元素求值 x np.arange(5) p(x) array([ 1., 6., 17., 34., 57.]) x np.arange(6).reshape(3,2) p(x) array([[ 1., 6.], [17., 34.], [57., 86.]])替换substitution把多项式代入x并展开若把多项式视为函数这就是函数复合 p(p) Polynomial([ 6., 16., 36., 36., 27.], domain[-1., 1.], window[-1., 1.], symbolx)求根 p.roots() array([-0.33333333-0.47140452j, -0.333333330.47140452j])roots()返回复根数组。从源码看它调用各类型的_roots后通过pu.mapdomain把根从window映射回domain见 _polybase.py#L900-L913文档字符串也提醒根越偏离domain区间精度越低。隐式类型转换与混合运算规则不必总是显式构造Polynomial实例元组、列表、数组和标量会在算术运算中被自动转换为系数参与计算 p [1, 2, 3] Polynomial([2., 4., 6.], domain[-1., 1.], window[-1., 1.], symbolx) [1, 2, 3] * p Polynomial([ 1., 4., 10., 12., 9.], domain[-1., 1.], window[-1., 1.], symbolx) p / 2 Polynomial([0.5, 1. , 1.5], domain[-1., 1.], window[-1., 1.], symbolx)domain、window 或类型不同的多项式不能直接混合运算 from numpy.polynomial import Chebyshev as T p P([1], domain[0,1]) TypeError: Domains differ p P([1], window[0,1]) TypeError: Windows differ p T([1]) TypeError: Polynomial types differ这些异常正是由_get_coefficients中的兼容性检查抛出的见 _polybase.py#L256-L290先检查类型是否一致再检查domain、window、symbol是否分别相同。利用替换实现类间转换不同类型不能直接做算术但可以用作替换substitution实现转换。事实上convert正是通过“把kind类的恒等多项式代入自身”来实现类型、domain、window 三重转换的见 convert 实现 p(T([0, 1])) Chebyshev([2.5, 2. , 1.5], domain[-1., 1.], window[-1., 1.], symbolx)这得到多项式p的 Chebyshev 形式。其原理是 $T_1(x) x$代入x不改变原多项式但所有乘除都在 Chebyshev 级数下进行因此结果类型为Chebyshev。不可变性与增强赋值所有多项式实例设计为不可变因此、-等增强赋值操作以及其他会破坏不可变性的功能都被有意未实现unimplemented。基类文档也声明__hash__ None并关闭了 ndarray 子类的 ufunc 与运算符交互见 _polybase.py#L70-L74。微积分微分与积分多项式实例支持积分与微分注意积分默认的下限与积分常数均为 0但都可自定义 from numpy.polynomial import Polynomial as P p P([2, 6]) p.integ() # 积分一次 Polynomial([0., 2., 3.], domain[-1., 1.], window[-1., 1.], symbolx) p.integ(2) # 积分两次 Polynomial([0., 0., 1., 1.], domain[-1., 1.], window[-1., 1.], symbolx) p.integ(lbnd-1) # 指定积分下限 Polynomial([-1., 2., 3.], domain[-1., 1.], window[-1., 1.], symbolx) p.integ(lbnd-1, k1) # 下限为 -1 且积分常数为 1 Polynomial([0., 2., 3.], domain[-1., 1.], window[-1., 1.], symbolx)m积分次数非负整数k积分常数数组k[0]应用于第一次积分、k[1]应用于第二次依此类推缺失值补零见 integ 签名。微分更简单唯一参数是微分次数 p P([1, 2, 3]) p.deriv(1) Polynomial([2., 6.], domain[-1., 1.], window[-1., 1.], symbolx) p.deriv(2) Polynomial([6.], domain[-1., 1.], window[-1., 1.], symbolx)值得一提的实现细节integ与deriv都先通过mapparms()获取domain→window的线性映射参数(off, scl)积分下限与微分结果都会按映射关系做缩放1./scl保证结果仍在该实例的domain上语义正确见 _polybase.py#L870-L898。其他构造方式fromroots、convert、basis 与 cast按系数构造只是获得多项式实例的一种方式还可以按根构造、从其他类型转换、以及通过最小二乘拟合构造拟合单独成节。按根构造与类型转换 from numpy.polynomial import Polynomial as P from numpy.polynomial import Chebyshev as T p P.fromroots([1, 2, 3]) p Polynomial([-6., 11., -6., 1.], domain[-1., 1.], window[-1., 1.], symbolx) p.convert(kindT) Chebyshev([-9. , 11.75, -3. , 0.25], domain[-1., 1.], window[-1., 1.], symbolx)fromroots构造的是 $\prod_i (x - r_i)$ 的展开其内部会先把根从domain映射到window再生成系数见 _polybase.py#L1037-L1077。convert方法还能同时转换 domain 与 window p.convert(kindT, domain[0, 1]) Chebyshev([-2.4375 , 2.96875, -0.5625 , 0.03125], domain[0., 1.], window[-1., 1.], symbolx) p.convert(kindP, domain[0, 1]) Polynomial([-1.875, 2.875, -1.125, 0.125], domain[0., 1.], window[-1., 1.], symbolx)basis 与 cast 类方法NumPy 1.7.0 起新增basis与cast类方法。cast与convert类似basis返回指定次数的基多项式 P.basis(3) Polynomial([0., 0., 0., 1.], domain[-1., 1.], window[-1., 1.], symbolx) T.cast(p) Chebyshev([-9. , 11.75, -3. , 0.25], domain[-1., 1.], window[-1., 1.], symbolx)basis的实现是把deg位置置 1、其余置 0 的系数交给构造器见 _polybase.py#L1114-L1151cast则委托给实例方法convert见 _polybase.py#L1153-L1191。类型转换的精度警告类型之间的转换虽然有用但不建议常规使用把 50 次的 Chebyshev 级数转换为同次幂级数时数值精度损失可能使求值结果近乎随机。基类的convert文档也明确提示“域和类型之间的转换可能产生数值上病态的级数”见 _polybase.py#L802-L806。拟合Fittingdomain 与 window 的用武之地拟合正是domain与window两个属性存在的根本原因。问题背景如下$[-1, 1]$ 区间内前几阶 Chebyshev 多项式是位于 ±1 之间的等波纹equiripple函数形态良好但同样的曲线放到 $[-2, 2]$ 区间良好区域就缩小到几乎可忽略。机制用 Chebyshev 多项式拟合时我们希望利用 $x$ 位于 $[-1, 1]$ 的“良好区域”——这正是window指定的区间。但待拟合数据几乎不可能恰好全部落在这个区间内于是用domain指定数据点所在的区间。拟合时先将domain通过线性变换映射到window再对映射后的数据点做常规最小二乘拟合。拟合结果的window与domain会随返回的级数一起保存并在后续求值、微分等操作中自动生效。若调用时未指定拟合例程会使用默认window和覆盖全部数据点的最小domain。线性映射的具体实现映射的数学核心在 polyutils.py 的三个函数中getdomain(x)返回覆盖数据的最小区间实数据为[min, max]复数据为包含全部点、边与坐标轴对齐的最小矩形的左下/右上角见 polyutils.py#L194-L239mapparms(old, new)返回线性映射L(x) off scl*x的参数满足old[i] → new[i]见 polyutils.py#L241-L286mapdomain(x, old, new)把映射应用到数据点x见 polyutils.py#L288-L323。fit的流程在 _polybase.py#L1014-L1034 中清晰可见domain缺省时由pu.getdomain(x)自动推导若数据退化为单点则自动扩展 ±1然后xnew pu.mapdomain(x, domain, window)映射数据点再交给具体类的_fit完成最小二乘。对含噪正弦数据的拟合示例 import numpy as np import matplotlib.pyplot as plt from numpy.polynomial import Chebyshev as T np.random.seed(11) x np.linspace(0, 2*np.pi, 20) y np.sin(x) np.random.normal(scale.1, sizex.shape) p T.fit(x, y, 5) plt.plot(x, y, o) xx, yy p.linspace() plt.plot(xx, yy, lw2) p.domain array([0. , 6.28318531]) p.window array([-1., 1.]) plt.show()注意因为调用时未显式指定domain拟合出的p.domain自动取为覆盖数据的最小区间[0, 2π]而p.window保持类默认的[-1, 1]——两者正是通过getdomain与mapparms联系起来的。fit 的完整参数fit类方法的签名见 _polybase.py#L946-L1012比文档示例更丰富x, y形如(M,)的样本点坐标deg拟合多项式的次数int表示包含直到deg次的所有项NumPy ≥ 1.11 还支持传入整数列表来指定要包含的特定次数项domainNone时自动选覆盖x的最小domain[]时用类默认domain[]选项 1.5.0 新增rcond相对条件数小于该值相对最大奇异值的奇异值被忽略默认len(x)*eps约 2e-16 量级fullTrue时额外返回[resid, rank, sv, rcond]诊断信息残差平方和、缩放 Vandermonde 矩阵的数值秩、奇异值、rcond详见linalg.lstsqw权重数组w[i]作用于x[i]处的未平方残差y[i] - y_hat[i]逆方差加权时取w[i] 1/sigma(y[i])window返回级数的window默认类默认windowsymbol自变量符号默认x。若需要未缩放、未平移的基多项式系数可执行new_series.convert().coef。绘图的得力助手linspace拟合示例中的p.linspace()是linspace(n100, domainNone)方法返回domain默认用实例自身的domain上n个等距点的(x, y)对主要用作绘图辅助见 _polybase.py#L915-L943。实用补充方法除本文已展开的方法外ABCPolyBase还提供一批常用工具方法见基类degree()返回次数系数长度减 1不检查非零尾项配合trim()可得到真实次数trim(tol0)移除绝对值不超过tol的尾随系数返回新实例cutdeg(deg)/truncate(size)按次数/长度截断级数常用于最小二乘中丢弃高阶小系数copy()返回副本mapparms()返回当前domain→window的线性映射参数(off, scl)has_samecoef/has_samedomain/has_samewindow/has_sametype一致性检查对应算术运算前的兼容性判定identity()恒等多项式p(x) x也是convert的实现基础。这些能力都有对应的测试覆盖例如 numpy/polynomial/tests/test_classes.py 与 test_chebyshev.py 等按类型划分的测试文件可作为进一步研读各类型底层_fit、_roots等实现行为的入口。小结NumPy 多项式便捷类以ABCPolyBase抽象基类统一了六种正交多项式级数的对象模型系数按次数递增存放、实例不可变、完整支持 - * // % divmod **运算与求值/复合/求根并提供integ/deriv微积分、fromroots/convert/basis/cast多种构造路径。真正发挥这些类威力的是domain/window映射机制——它把任意区间上的数据线性映射到基函数表现良好的 $[-1,1]$ 区间再做最小二乘拟合并在返回的实例中持久保存映射关系使后续求值、微分自动一致。理解这一机制是使用Chebyshev.fit等高阶拟合能力时避免数值病态的关键。【免费下载链接】numpyThe fundamental package for scientific computing with Python.项目地址: https://gitcode.com/gh_mirrors/nu/numpy创作声明：本文部分内容由AI辅助生成（AIGC），仅供参考