跳到主要内容

NumPy 多项式

多项式是数学和科学计算中的基础工具，广泛应用于曲线拟合、信号处理、数值积分等领域。NumPy 提供了 numpy.polynomial 模块，支持多种多项式类型的创建、运算和拟合。本章将详细介绍多项式的使用方法。

多项式基础​

什么是多项式？​

多项式是由变量和系数通过加、减、乘和正整数次幂运算构成的代数表达式。一般形式为：

P(x)=c0+c1x+c2x2+⋯+cnxnP(x) = c_0 + c_1 x + c_2 x^2 + \cdots + c_n x^n

其中 c0,c1,…,cnc_0, c_1, \ldots, c_n 是系数，nn 是多项式的次数。

多项式类型​

NumPy 支持多种正交多项式：

多项式类型说明适用场景
Polynomial普通多项式（幂级数）通用计算
Chebyshev切比雪夫多项式数值逼近、滤波
Legendre勒让德多项式物理问题、球谐函数
Laguerre拉盖尔多项式量子力学
Hermite埃尔米特多项式量子力学、概率论

创建多项式​

从系数创建​

import numpy as np
from numpy.polynomial import Polynomial

# 从系数创建多项式
# 系数按升序排列：c0, c1, c2, ...
# P(x) = 1 + 2x + 3x^2
p = Polynomial([1, 2, 3])

print(f"多项式: {p}")
print(f"系数: {p.coef}")
print(f"定义域: {p.domain}")
print(f"窗口: {p.window}")

# 计算多项式在 x=2 处的值
x = 2
result = p(x)
print(f"\nP({x}) = {result}")  # 1 + 2*2 + 3*4 = 17

# 在多个点计算
x_values = np.array([0, 1, 2, 3])
print(f"P({x_values}) = {p(x_values)}")

理解 domain 和 window 属性​

NumPy 多项式有两个重要的属性：domain 和 window。这两个属性在多项式拟合中起着关键作用，但常常被忽视。

domain（定义域）：指定多项式适用的 x 值范围。当使用 fit 方法时，domain 会自动设置为数据的 x 值范围。

window（窗口）：定义多项式系数"标准化"后的 x 值范围，默认为 [-1, 1]。

为什么要引入这两个属性？这涉及到数值计算的稳定性问题。

当进行多项式拟合时，如果原始数据的 x 值范围很大（比如 [1000, 2000]），直接计算 xnx^n 会导致数值溢出或精度丢失。NumPy 的解决方案是：

  1. 将原始 x 值从 [domain[0], domain[1]] 线性映射到 [window[0], window[1]]
  2. 在标准化后的窗口内计算多项式系数
  3. 求值时自动进行逆变换
import numpy as np
from numpy.polynomial import Polynomial

# 演示 domain 和 window 的作用
x_data = np.array([0, 10, 20, 30, 40])  # 数据范围较大
y_data = np.array([1, 25, 81, 169, 289])  # y ≈ x^2 + 2x + 1

# 拟合多项式
p = Polynomial.fit(x_data, y_data, deg=2)

print(f"拟合后的 domain: {p.domain}")  # 自动设置为 [0, 40]
print(f"拟合后的 window: {p.window}")  # 默认 [-1, 1]

# 直接查看系数（在标准化窗口内的系数）
print(f"\n标准化系数: {p.coef}")

# 转换到原始数据域的系数
p_converted = p.convert(domain=[-1, 1], window=[-1, 1])
print(f"原始域系数: {p_converted.coef}")

# 验证：两种方式求值结果相同
test_x = 25
print(f"\n验证 P({test_x}):")
print(f"  原始多项式: {p(test_x):.2f}")
print(f"  转换后多项式: {p_converted(test_x):.2f}")

关键理解：fit 方法返回的多项式，其系数是基于标准化窗口计算的。如果需要得到"真实"的系数（对应原始数据域），需要使用 convert() 方法。

# 实际案例：温度传感器校准
temp_measured = np.array([20, 40, 60, 80, 100])  # 摄氏度
voltage = np.array([0.8, 1.6, 2.4, 3.2, 4.0])    # 伏特

# 拟合：电压 = f(温度)
calibration = Polynomial.fit(temp_measured, voltage, deg=1)

print("温度传感器校准:")
print(f"  domain: {calibration.domain}")  # [20, 100]

# 使用 convert 获取实际斜率和截距
cal_std = calibration.convert(domain=[-1, 1], window=[-1, 1])
print(f"  斜率 (V/°C): {cal_std.coef[1] / 40}")  # 注意域映射
print(f"  截距 (V): {cal_std.coef[0] - cal_std.coef[1] * 60 / 40}")

# 更直接的方式：使用 linspace 获取"真实"系数
x_test = np.array([0, 1])
y_test = calibration(x_test)
slope = y_test[1] - y_test[0]
intercept = y_test[0]
print(f"\n直接计算:")
print(f"  斜率 (V/°C): {slope:.4f}")
print(f"  截距 (V): {intercept:.4f}")

从根创建​

import numpy as np
from numpy.polynomial import Polynomial

# 从根创建多项式
# 根为 1, 2, 3 的多项式是 (x-1)(x-2)(x-3)
p = Polynomial.fromroots([1, 2, 3])

print(f"多项式: {p}")
print(f"展开形式: {p.coef}")

# 验证根
print(f"\n验证 P(1) = {p(1):.2f}")  # 接近 0
print(f"验证 P(2) = {p(2):.2f}")  # 接近 0
print(f"验证 P(3) = {p(3):.2f}")  # 接近 0

特殊多项式​

import numpy as np
from numpy.polynomial import Polynomial

# 零多项式
p_zero = Polynomial.zero()
print(f"零多项式: {p_zero}")

# 常数多项式
p_const = Polynomial.constant(5)
print(f"常数多项式: {p_const}")

# 恒等多项式 P(x) = x
p_identity = Polynomial.identity()
print(f"恒等多项式: {p_identity}")
print(f"P(10) = {p_identity(10)}")

多项式运算​

基本算术运算​

import numpy as np
from numpy.polynomial import Polynomial

p1 = Polynomial([1, 2, 3])  # 1 + 2x + 3x^2
p2 = Polynomial([4, 5])     # 4 + 5x

# 加法
add = p1 + p2
print(f"加法: {p1} + {p2} = {add}")

# 减法
sub = p1 - p2
print(f"减法: {p1} - {p2} = {sub}")

# 乘法
mul = p1 * p2
print(f"乘法: {p1} * {p2} = {mul}")

# 除法（返回商和余数）
quotient, remainder = p1 // p2
print(f"除法: {p1} // {p2} = {quotient}，余数 = {remainder}")

# 幂运算
pow_result = p1 ** 2
print(f"平方: ({p1})^2 = {pow_result}")

多项式复合​

import numpy as np
from numpy.polynomial import Polynomial

p1 = Polynomial([0, 1])    # P1(x) = x
p2 = Polynomial([1, 2, 1]) # P2(x) = 1 + 2x + x^2 = (1+x)^2

# 复合：P1(P2(x))
composed = p1(p2)
print(f"P1(P2(x)) = {composed}")

# 复合：P2(P1(x))
composed2 = p2(p1)
print(f"P2(P1(x)) = {composed2}")

导数和积分​

import numpy as np
from numpy.polynomial import Polynomial

p = Polynomial([1, 2, 3, 4])  # 1 + 2x + 3x^2 + 4x^3

# 导数
p_deriv = p.deriv()
print(f"原多项式: {p}")
print(f"一阶导数: {p_deriv}")

# 二阶导数
p_deriv2 = p.deriv(m=2)
print(f"二阶导数: {p_deriv2}")

# 积分
p_integ = p.integ()
print(f"\n积分（默认常数项为0）: {p_integ}")

# 积分并指定常数项
p_integ_const = p.integ(k=5)
print(f"积分（常数项=5）: {p_integ_const}")

# 验证：导数和积分互为逆运算
print(f"\n验证: 积分后再求导 = {p_integ.deriv()}")

多项式求值与求根​

求值​

import numpy as np
from numpy.polynomial import Polynomial

p = Polynomial([1, 2, 3])  # 1 + 2x + 3x^2

# 单点求值
print(f"P(2) = {p(2)}")

# 多点求值
x = np.linspace(-2, 2, 5)
print(f"P({x}) = {p(x)}")

# 使用 __call__ 方法（等价于直接调用）
values = p.__call__(x)
print(f"使用 __call__: {values}")

求根​

import numpy as np
from numpy.polynomial import Polynomial

# 创建多项式并求根
p = Polynomial([6, -5, 1])  # 6 - 5x + x^2 = (x-2)(x-3)
roots = p.roots()

print(f"多项式: {p}")
print(f"根: {roots}")

# 验证
for root in roots:
    print(f"P({root:.2f}) = {p(root):.2e}")

从值求系数​

import numpy as np
from numpy.polynomial import Polynomial

# 已知 x 和 y 值，求多项式系数
x = np.array([0, 1, 2, 3])
y = np.array([1, 6, 17, 34])  # y = 1 + 2x + 3x^2

# 使用 polyfit 拟合
p = Polynomial.fit(x, y, deg=2)

print(f"拟合多项式: {p}")
print(f"拟合系数: {p.convert().coef}")  # 转换到标准形式

# 验证
print(f"\n验证:")
for xi, yi in zip(x, y):
    print(f"  P({xi}) = {p(xi):.2f}, 期望值 = {yi}")

多项式拟合​

多项式拟合是数据分析中最常用的技术之一。理解拟合的本质、如何选择合适的次数、以及如何避免常见陷阱，对于实际应用至关重要。

拟合的数学原理​

给定 mm 个数据点 (xi,yi)(x_i, y_i)，寻找 nn 次多项式 P(x)=∑k=0nckxkP(x) = \sum_{k=0}^{n} c_k x^k，使得残差平方和最小：

min⁡c0,c1,…,cn∑i=1m(yi−P(xi))2\min_{c_0, c_1, \ldots, c_n} \sum_{i=1}^{m} (y_i - P(x_i))^2

这是一个线性最小二乘问题。令 VV 为 Vandermonde 矩阵：

V=[1x1x12⋯x1n1x2x22⋯x2n⋮⋮⋮⋱⋮1xmxm2⋯xmn]V = \begin{bmatrix} 1 & x_1 & x_1^2 & \cdots & x_1^n \\ 1 & x_2 & x_2^2 & \cdots & x_2^n \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & x_m & x_m^2 & \cdots & x_m^n \end{bmatrix}

则系数向量 c=(c0,c1,…,cn)T\mathbf{c} = (c_0, c_1, \ldots, c_n)^T 的最优解为：

c=(VTV)−1VTy\mathbf{c} = (V^T V)^{-1} V^T \mathbf{y}

NumPy 内部使用更稳定的数值方法（如 QR 分解或 SVD）来求解。

基本拟合​

import numpy as np
from numpy.polynomial import Polynomial

# 生成带噪声的数据
np.random.seed(42)
x = np.linspace(-3, 3, 30)
y_true = 1 + 2*x - x**2  # 真实函数
y = y_true + np.random.normal(0, 0.5, len(x))  # 添加噪声

# 多项式拟合
p_fit = Polynomial.fit(x, y, deg=2)

print(f"拟合多项式: {p_fit}")

# 预测
x_pred = np.linspace(-3, 3, 100)
y_pred = p_fit(x_pred)

# 计算拟合误差
residual = np.sum((y - p_fit(x))**2)
print(f"残差平方和: {residual:.4f}")

过拟合与欠拟合​

选择多项式次数是拟合中的关键决策。次数太低导致欠拟合，次数太高导致过拟合。

import numpy as np
from numpy.polynomial import Polynomial

np.random.seed(42)

# 生成数据：真实关系是二次函数
x_train = np.linspace(-2, 2, 20)
y_true = 1 + 2*x - 0.5*x**2
y_train = y_true[:len(x_train)] + 0.3 * np.random.randn(len(x_train))

# 测试数据
x_test = np.linspace(-2.5, 2.5, 30)  # 包含外推区域
y_test_true = 1 + 2*x_test - 0.5*x_test**2
y_test = y_test_true + 0.3 * np.random.randn(len(x_test))

print("拟合分析：")
print("-" * 65)
print(f"{'次数':>4} | {'训练误差':>12} | {'测试误差':>12} | {'评价':>15}")
print("-" * 65)

for deg in [1, 2, 4, 8, 15]:
    p = Polynomial.fit(x_train, y_train, deg=deg)
    
    train_error = np.mean((y_train - p(x_train))**2)
    test_error = np.mean((y_test - p(x_test))**2)
    
    if deg == 1:
        status = "欠拟合"
    elif deg == 2:
        status = "最佳"
    elif test_error > train_error * 10:
        status = "严重过拟合"
    else:
        status = "轻度过拟合"
    
    print(f"{deg:4d} | {train_error:12.6f} | {test_error:12.6f} | {status:>15}")

判断过拟合的方法：

  1. 训练误差远小于测试误差：典型过拟合标志
  2. 系数绝对值很大：多项式试图用大的系数"抵消"相互影响
  3. 曲线剧烈振荡：高次多项式在数据点间震荡

使用交叉验证选择最佳次数​

import numpy as np
from numpy.polynomial import Polynomial

def cross_validate_degree(x, y, max_deg=10, k_folds=5):
    """使用 k 折交叉验证选择最佳多项式次数"""
    n = len(x)
    fold_size = n // k_folds
    
    best_deg = 1
    best_cv_error = float('inf')
    
    print("交叉验证结果:")
    for deg in range(1, max_deg + 1):
        cv_errors = []
        
        for fold in range(k_folds):
            # 划分训练集和验证集
            val_start = fold * fold_size
            val_end = val_start + fold_size
            
            x_val = x[val_start:val_end]
            y_val = y[val_start:val_end]
            x_train = np.concatenate([x[:val_start], x[val_end:]])
            y_train = np.concatenate([y[:val_start], y[val_end:]])
            
            # 拟合并计算验证误差
            p = Polynomial.fit(x_train, y_train, deg=deg)
            error = np.mean((y_val - p(x_val))**2)
            cv_errors.append(error)
        
        mean_cv_error = np.mean(cv_errors)
        print(f"  次数 {deg}: CV 误差 = {mean_cv_error:.6f}")
        
        if mean_cv_error < best_cv_error:
            best_cv_error = mean_cv_error
            best_deg = deg
    
    return best_deg

# 示例
np.random.seed(42)
x = np.linspace(-3, 3, 50)
y = 1 + 2*x - x**2 + 0.3*np.random.randn(len(x))

best_deg = cross_validate_degree(x, y, max_deg=10)
print(f"\n推荐的多项式次数: {best_deg}")

正则化：岭回归拟合​

对于高次多项式，可以使用正则化来防止过拟合：

import numpy as np

def ridge_polyfit(x, y, deg, alpha=0.1):
    """
    使用岭回归进行多项式拟合
    
    参数:
        x, y: 数据点
        deg: 多项式次数
        alpha: 正则化参数（越大，系数越平滑）
    """
    from numpy.polynomial import polynomial as P
    
    # 构建 Vandermonde 矩阵
    V = P.polyvander(x, deg)
    
    # 岭回归: (V^T V + alpha*I)^{-1} V^T y
    I = np.eye(deg + 1)
    I[0, 0] = 0  # 通常不对常数项正则化
    
    coef = np.linalg.solve(V.T @ V + alpha * I, V.T @ y)
    return coef

# 示例：过拟合情况下使用正则化
np.random.seed(42)
x = np.linspace(-2, 2, 15)
y = 1 + 0.5*x + 0.1*np.random.randn(len(x))

# 普通拟合（过拟合）
from numpy.polynomial import Polynomial
p_normal = Polynomial.fit(x, y, deg=10)

# 岭回归拟合
coef_ridge = ridge_polyfit(x, y, deg=10, alpha=0.5)
p_ridge = Polynomial(coef_ridge)

# 比较系数大小
print("系数绝对值对比:")
print(f"  普通拟合最大系数: {np.max(np.abs(p_normal.convert().coef)):.2f}")
print(f"  岭回归最大系数: {np.max(np.abs(coef_ridge)):.2f}")

# 外推测试
x_test = 2.5
print(f"\n外推到 x={x_test}:")
print(f"  普通拟合: {p_normal(x_test):.2f}")
print(f"  岭回归: {p_ridge(x_test):.2f}")
# 岭回归的外推更稳定

不同次数的拟合比较​

import numpy as np
from numpy.polynomial import Polynomial

# 生成复杂数据
np.random.seed(42)
x = np.linspace(-2, 2, 50)
y = np.sin(2*x) + 0.1 * np.random.randn(len(x))

print("不同次数多项式拟合比较:")
for deg in [1, 3, 5, 7]:
    p = Polynomial.fit(x, y, deg=deg)
    residual = np.sum((y - p(x))**2)
    print(f"  次数 {deg}: 残差平方和 = {residual:.4f}")

使用 Vandermonde 矩阵​

import numpy as np
from numpy.polynomial import polynomial as P

# 生成 Vandermonde 矩阵
x = np.array([1, 2, 3, 4])
deg = 3

# Vandermonde 矩阵：每一列是 x 的不同次幂
V = P.polyvander(x, deg)
print(f"Vandermonde 矩阵:\n{V}")

# 用最小二乘法求解系数
y = np.array([1, 4, 9, 16])  # y = x^2
coef = np.linalg.lstsq(V, y, rcond=None)[0]
print(f"\n拟合系数: {coef}")

切比雪夫多项式​

切比雪夫多项式在数值分析中非常重要，具有良好的数值稳定性。理解切比雪夫多项式对于高精度数值计算至关重要。

为什么选择切比雪夫多项式？​

普通多项式（幂级数）在高次拟合时存在严重的数值问题。考虑拟合区间 [-1, 1] 上的函数，使用普通多项式 P(x)=∑k=0nckxkP(x) = \sum_{k=0}^{n} c_k x^k：

  • 当 xx 接近 ±1\pm 1 时，xnx^n 的值在区间内变化很大
  • 高次项系数往往很小，但乘以 xnx^n 后可能很大
  • 基函数 1,x,x2,x3,…1, x, x^2, x^3, \ldots 之间高度相关，导致病态的线性方程组

切比雪夫多项式解决了这个问题。切比雪夫多项式 Tn(x)T_n(x) 定义为：

Tn(x)=cos⁡(narccos⁡(x)),x∈[−1,1]T_n(x) = \cos(n \arccos(x)), \quad x \in [-1, 1]

前几项切比雪夫多项式为：

  • T0(x)=1T_0(x) = 1
  • T1(x)=xT_1(x) = x
  • T2(x)=2x2−1T_2(x) = 2x^2 - 1
  • T3(x)=4x3−3xT_3(x) = 4x^3 - 3x
  • Tn(x)=2x⋅Tn−1(x)−Tn−2(x)T_n(x) = 2x \cdot T_{n-1}(x) - T_{n-2}(x)

切比雪夫多项式的主要优势：

1. 正交性：在区间 [−1,1][-1, 1] 上，切比雪夫多项式关于权重函数 11−x2\frac{1}{\sqrt{1-x^2}} 正交：

∫−11Tm(x)Tn(x)1−x2dx={0m≠nπm=n=0π2m=n≠0\int_{-1}^{1} \frac{T_m(x) T_n(x)}{\sqrt{1-x^2}} dx = \begin{cases} 0 & m \neq n \\ \pi & m = n = 0 \\ \frac{\pi}{2} & m = n \neq 0 \end{cases}

2. 等振荡性质：切比雪夫多项式在 [−1,1][-1, 1] 上的极值点值相等，都为 ±1\pm 1。这意味着每一项对结果的影响更均匀。

3. 最佳逼近：对于给定的连续函数，切比雪夫展开往往比幂级数展开收敛更快。

创建切比雪夫多项式​

import numpy as np
from numpy.polynomial import Chebyshev

# 从系数创建
c = Chebyshev([1, 2, 3])  # T0(x) + 2*T1(x) + 3*T2(x)
print(f"切比雪夫多项式: {c}")

# 计算值
x = 0.5
print(f"C({x}) = {c(x)}")

# 与普通多项式转换
p = c.convert(kind=np.polynomial.Polynomial)
print(f"转换为普通多项式: {p}")

切比雪夫逼近的收敛性​

import numpy as np
from numpy.polynomial import Chebyshev, Polynomial

# 比较切比雪夫和普通多项式逼近
x = np.linspace(-1, 1, 200)
y = np.exp(x)  # 要逼近的函数

print("逼近 e^x 的误差比较:")
print("次数 | 普通多项式误差 | 切比雪夫误差")
print("-" * 45)

for deg in [2, 4, 6, 8, 10]:
    # 普通多项式
    p_fit = Polynomial.fit(x, y, deg=deg)
    p_error = np.max(np.abs(y - p_fit(x)))
    
    # 切比雪夫多项式
    c_fit = Chebyshev.fit(x, y, deg=deg)
    c_error = np.max(np.abs(y - c_fit(x)))
    
    print(f"  {deg:2d} |     {p_error:.2e}    |   {c_error:.2e}")

# 可以看到：对于同样的次数，切比雪夫多项式的逼近误差通常更小

Runge 现象与切比雪夫节点​

等距插值点会导致 Runge 现象——在区间边界处出现剧烈振荡。切比雪夫节点可以避免这个问题：

import numpy as np
from numpy.polynomial import Polynomial

def runge_function(x):
    """Runge 函数：等距节点插值的经典反例"""
    return 1 / (1 + 25 * x**2)

# 等距节点 vs 切比雪夫节点
n = 15

# 等距节点
x_equidistant = np.linspace(-1, 1, n)
y_equi = runge_function(x_equidistant)
p_equi = Polynomial.fit(x_equidistant, y_equi, deg=n-1)

# 切比雪夫节点
k = np.arange(1, n + 1)
x_chebyshev = np.cos((2 * k - 1) * np.pi / (2 * n))  # 切比雪夫节点公式
y_cheb = runge_function(x_chebyshev)
p_cheb = Polynomial.fit(x_chebyshev, y_cheb, deg=n-1)

# 评估
x_test = np.linspace(-1, 1, 1000)
y_true = runge_function(x_test)

error_equi = np.max(np.abs(y_true - p_equi(x_test)))
error_cheb = np.max(np.abs(y_true - p_cheb(x_test)))

print(f"等距节点插值最大误差: {error_equi:.4f}")
print(f"切比雪夫节点插值最大误差: {error_cheb:.4f}")
# 切比雪夫节点的误差显著更小，没有边界振荡问题

切比雪夫逼近​

import numpy as np
from numpy.polynomial import Chebyshev

# 使用切比雪夫多项式逼近函数
x = np.linspace(-1, 1, 100)
y = np.exp(x)  # 要逼近的函数

# 用不同次数的切比雪夫多项式逼近
for deg in [2, 4, 6, 8]:
    c_fit = Chebyshev.fit(x, y, deg=deg)
    max_error = np.max(np.abs(y - c_fit(x)))
    print(f"次数 {deg}: 最大误差 = {max_error:.6f}")

切比雪夫根（最佳插值点）​

import numpy as np
from numpy.polynomial import Chebyshev

# 切比雪夫节点：在区间上均匀分布的根
n = 10
# 在 [-1, 1] 区间
cheb_roots = Chebyshev.basis(n).roots()
print(f"切比雪夫根: {cheb_roots}")

# 转换到 [a, b] 区间
a, b = 0, 10
cheb_roots_ab = (b - a) / 2 * cheb_roots + (a + b) / 2
print(f"在 [{a}, {b}] 区间的根: {cheb_roots_ab}")

其他正交多项式​

勒让德多项式​

import numpy as np
from numpy.polynomial import Legendre

# 创建勒让德多项式
l = Legendre([1, 2, 3])  # P0(x) + 2*P1(x) + 3*P2(x)
print(f"勒让德多项式: {l}")

# 勒让德多项式在 [-1, 1] 上正交
x = np.linspace(-1, 1, 100)
P0 = Legendre.basis(0)(x)  # 常数 1
P1 = Legendre.basis(1)(x)  # x
P2 = Legendre.basis(2)(x)  # (3x^2 - 1)/2

# 验证正交性（积分应该接近 0）
integral_01 = np.trapezoid(P0 * P1, x)
integral_02 = np.trapezoid(P0 * P2, x)
integral_12 = np.trapezoid(P1 * P2, x)
print(f"\n正交性验证:")
print(f"∫P0*P1 dx = {integral_01:.6f}")
print(f"∫P0*P2 dx = {integral_02:.6f}")
print(f"∫P1*P2 dx = {integral_12:.6f}")

拉盖尔和埃尔米特多项式​

import numpy as np
from numpy.polynomial import Laguerre, Hermite

# 拉盖尔多项式（在 [0, ∞) 上正交，权重函数 e^(-x)）
lag = Laguerre([1, 1, 1])
print(f"拉盖尔多项式: {lag}")

# 埃尔米特多项式（在 (-∞, ∞) 上正交，权重函数 e^(-x^2)）
herm = Hermite([1, 0, 1])
print(f"埃尔米特多项式: {herm}")

实际应用​

应用1：曲线拟合与模型选择​

在实际数据分析中，多项式拟合常用于建立变量间的关系模型。关键是要选择合适的模型复杂度。

import numpy as np
from numpy.polynomial import Polynomial

# 案例：实验数据建模
# 假设测量了一组温度与电阻的关系
temperature = np.array([0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100])
resistance = np.array([100, 103.9, 107.7, 111.5, 115.3, 119.0, 122.6, 126.2, 129.7, 133.2, 136.6])

# 尝试不同次数的拟合
print("电阻-温度关系拟合分析:")
print("-" * 50)

best_deg = 1
best_aic = float('inf')

for deg in range(1, 5):
    p = Polynomial.fit(temperature, resistance, deg=deg)
    
    # 计算残差
    residual = resistance - p(temperature)
    rss = np.sum(residual**2)  # 残差平方和
    n = len(temperature)
    
    # 使用 AIC 准则选择模型（平衡拟合优度和模型复杂度）
    # AIC = n * ln(RSS/n) + 2*k，k = deg + 1
    k = deg + 1
    aic = n * np.log(rss / n) + 2 * k
    
    print(f"次数 {deg}: RSS = {rss:.2f}, AIC = {aic:.2f}")
    
    if aic < best_aic:
        best_aic = aic
        best_deg = deg

print(f"\n根据 AIC 准则，推荐次数: {best_deg}")

# 使用最佳模型
p_best = Polynomial.fit(temperature, resistance, deg=best_deg)
p_std = p_best.convert(domain=[-1, 1], window=[-1, 1])

print(f"\n拟合多项式系数: {p_std.coef}")

# 预测新温度下的电阻
temp_pred = np.array([25, 55, 85])
resist_pred = p_best(temp_pred)
print(f"\n预测:")
for t, r in zip(temp_pred, resist_pred):
    print(f"  {t}°C: {r:.1f} Ω")

应用2：数值积分​

多项式的一个重要应用是数值积分。通过多项式拟合函数，可以利用多项式易积分的性质计算复杂函数的积分。

import numpy as np
from numpy.polynomial import Polynomial, Chebyshev

def poly_integrate(func, a, b, n_points=50, deg=10, use_chebyshev=False):
    """
    使用多项式逼近进行数值积分
    
    参数:
        func: 被积函数
        a, b: 积分区间
        n_points: 采样点数
        deg: 多项式次数
        use_chebyshev: 是否使用切比雪夫多项式
    """
    if use_chebyshev:
        # 切比雪夫节点（更稳定）
        k = np.arange(1, n_points + 1)
        x = (b - a) / 2 * np.cos((2 * k - 1) * np.pi / (2 * n_points)) + (a + b) / 2
        x = np.sort(x)  # 排序以便积分计算
        y = func(x)
        p = Chebyshev.fit(x, y, deg=deg)
    else:
        # 等距节点
        x = np.linspace(a, b, n_points)
        y = func(x)
        p = Polynomial.fit(x, y, deg=deg)
    
    # 多项式积分
    p_integ = p.integ()
    
    # 定积分
    result = p_integ(b) - p_integ(a)
    return result

# 测试：积分 sin(x) 在 [0, π]
result = poly_integrate(np.sin, 0, np.pi)
true_value = 2.0
print(f"积分 sin(x) 在 [0, π]:")
print(f"  多项式法: {result:.10f}")
print(f"  真实值: {true_value:.10f}")
print(f"  误差: {abs(result - true_value):.2e}")

# 使用切比雪夫多项式
result_cheb = poly_integrate(np.sin, 0, np.pi, use_chebyshev=True)
print(f"\n切比雪夫法: {result_cheb:.10f}")
print(f"切比雪夫误差: {abs(result_cheb - true_value):.2e}")

应用3：牛顿法求根​

多项式求根在实际问题中广泛应用，如求解方程、寻找最优值等。

import numpy as np
from numpy.polynomial import Polynomial

def newton_find_root(p, x0, tol=1e-10, max_iter=100):
    """
    使用牛顿法求多项式的根
    
    牛顿法迭代公式: x_{n+1} = x_n - P(x_n) / P'(x_n)
    """
    p_deriv = p.deriv()
    
    x = x0
    for i in range(max_iter):
        fx = p(x)
        fpx = p_deriv(x)
        
        if abs(fpx) < 1e-15:
            print("警告：导数接近零，可能无法继续")
            break
        
        x_new = x - fx / fpx
        
        if abs(x_new - x) < tol:
            return x_new, i + 1
        
        x = x_new
    
    return x, max_iter

# 示例：求 x^3 - 2x - 5 = 0 的根
p = Polynomial([-5, -2, 0, 1])  # x^3 - 2x - 5
print(f"方程: {p} = 0")

# 使用 NumPy 内置方法
roots = p.roots()
print(f"\nNumPy 求根结果: {roots}")

# 使用牛顿法（需要好的初值）
initial_guess = 2.0
root_newton, iterations = newton_find_root(p, initial_guess)
print(f"\n牛顿法求根:")
print(f"  初值: {initial_guess}")
print(f"  根: {root_newton:.10f}")
print(f"  迭代次数: {iterations}")
print(f"  验证 P({root_newton:.6f}) = {p(root_newton):.2e}")

应用4：信号平滑与滤波​

多项式局部拟合是一种有效的信号平滑方法，称为 Savitzky-Golay 滤波器的思想基础。

import numpy as np
from numpy.polynomial import Polynomial

def savitzky_golay_simple(signal, window_size=11, poly_deg=2):
    """
    简化的 Savitzky-Golay 滤波器
    
    原理：在每个窗口内用多项式拟合，取多项式在中心点的值作为平滑结果
    """
    half_window = window_size // 2
    smoothed = np.zeros_like(signal, dtype=float)
    n = len(signal)
    
    for i in range(n):
        # 确定窗口边界
        left = max(0, i - half_window)
        right = min(n, i + half_window + 1)
        
        # 调整窗口使其对称（边界处理）
        actual_left = i - (right - i - 1)
        actual_right = i + (i - left) + 1
        left = max(0, actual_left)
        right = min(n, actual_right)
        
        # 窗口内拟合
        x_window = np.arange(left, right)
        y_window = signal[left:right]
        
        if len(x_window) > poly_deg:
            p = Polynomial.fit(x_window, y_window, deg=poly_deg)
            smoothed[i] = p(i)
        else:
            smoothed[i] = signal[i]
    
    return smoothed

# 示例：平滑含噪声信号
np.random.seed(42)
t = np.linspace(0, 4*np.pi, 200)
signal_clean = np.sin(t)
signal_noisy = signal_clean + 0.3 * np.random.randn(len(t))

# 平滑处理
signal_smoothed = savitzky_golay_simple(signal_noisy, window_size=15, poly_deg=3)

# 评估效果
noise_before = np.std(signal_noisy - signal_clean)
noise_after = np.std(signal_smoothed - signal_clean)
print(f"平滑效果:")
print(f"  原始噪声标准差: {noise_before:.4f}")
print(f"  平滑后标准差: {noise_after:.4f}")
print(f"  噪声降低: {(1 - noise_after/noise_before)*100:.1f}%")

应用5：寻找函数极值​

import numpy as np
from numpy.polynomial import Polynomial

def find_extrema(p, domain=(-10, 10)):
    """
    找到多项式在给定区间内的所有极值点
    
    极值点是导数为零的点
    """
    # 导数
    p_deriv = p.deriv()
    
    # 导数的根即极值点
    roots = p_deriv.roots()
    
    # 过滤实根和在区间内的根
    extrema = []
    p_deriv2 = p_deriv.deriv()
    
    for root in roots:
        if np.isreal(root):
            x = root.real
            if domain[0] <= x <= domain[1]:
                # 判断极值类型
                second_deriv = p_deriv2(x)
                extrema_type = "极小值" if second_deriv > 0 else "极大值"
                extrema.append((x, p(x), extrema_type))
    
    return extrema

# 示例：分析多项式的极值
p = Polynomial([1, -3, 3, -1])  # 1 - 3x + 3x^2 - x^3
print(f"多项式: {p}")

extrema = find_extrema(p, domain=(-5, 5))
print(f"\n极值点分析:")
for x, y, etype in extrema:
    print(f"  x = {x:.4f}, y = {y:.4f}, 类型: {etype}")

# 找全局最值
x_eval = np.linspace(-10, 10, 1000)
y_eval = p(x_eval)
print(f"\n全局最值:")
print(f"  全局最大值: {np.max(y_eval):.4f} (在 x = {x_eval[np.argmax(y_eval)]:.4f})")
print(f"  全局最小值: {np.min(y_eval):.4f} (在 x = {x_eval[np.argmin(y_eval)]:.4f})")

多项式与旧 API 对比​

NumPy 有两套多项式 API：旧版 np.poly1d 和新版 np.polynomial：

import numpy as np

# 旧 API（不推荐，但仍可用）
p_old = np.poly1d([3, 2, 1])  # 系数按降序：3x^2 + 2x + 1
print(f"旧 API: {p_old}")

# 新 API（推荐）
from numpy.polynomial import Polynomial
p_new = Polynomial([1, 2, 3])  # 系数按升序：1 + 2x + 3x^2
print(f"新 API: {p_new}")

# 注意系数顺序相反！

新版 API 的优势：

  • 系数顺序更符合数学习惯（升序）
  • 支持多种正交多项式
  • 更好的数值稳定性
  • 更完整的功能（拟合、积分、导数等）

NumPy 2.0 多项式变化​

NumPy 2.0 对多项式的表示方法进行了更新：

import numpy as np
from numpy.polynomial import Polynomial

# NumPy 2.0 中，多项式的字符串表示包含 domain 信息
p = Polynomial([1, 2, 3], domain=[0, 1])

# 字符串表示现在更详细
print(f"NumPy 2.0 多项式表示: {p}")
# 输出格式: 1.0 + 2.0 x + 3.0 x^2
# 如果设置了自定义 domain，会显示转换公式

# 带自定义 domain 的多项式
p_domain = Polynomial([1, 2, 3], domain=[0.1, 0.2])
print(f"\n带 domain 的多项式: {p_domain}")
# NumPy 2.0 会显示 domain 信息

NumPy 2.0 主要变化：

  1. 表示方法更新：多项式的 str 和 repr 表示现在包含 domain 信息
  2. 一致的输出：纯文本和 LaTeX 表示现在保持一致
  3. 更好的类型提示：改进了类型注解支持
# 示例：NumPy 2.0 中 Polynomial 的完整表示
p = Polynomial([1, 1], domain=[0.1, 0.2])
print(f"完整表示: {p}")
# 显示: 1.0 + 1.0 (-3.0000000000000004 + 20.0 x)
# 这表示在 domain [0.1, 0.2] 上的变换

函数参考​

Polynomial 类方法​

方法说明
__call__(x)求值
roots()求根
deriv(m=1)导数
integ(m=1, k=0)积分
convert(domain, window, kind)类型转换
fit(x, y, deg)数据拟合
fromroots(roots)从根创建
trunc(size)截断到指定次数
copy()复制

工厂方法​

方法说明
Polynomial.basis(deg)创建指定次数的基多项式
Polynomial.identity()恒等多项式 P(x) = x
Polynomial.zero()零多项式
Polynomial.constant(c)常数多项式

小结​

本章介绍了 NumPy 多项式模块的核心功能：

  1. 多项式创建：从系数、从根、特殊多项式
  2. 基本运算：加减乘除、幂、复合
  3. 微积分：导数、积分
  4. 求值与求根：在指定点求值、求解方程
  5. 数据拟合：多项式回归、切比雪夫逼近
  6. 正交多项式：切比雪夫、勒让德、拉盖尔、埃尔米特
  7. 实际应用：曲线拟合、数值积分、信号处理

多项式是数值计算的基础工具，理解其使用方法对于科学计算和数据分析至关重要。

练习​

  1. 创建一个多项式，找到它在区间 [-2, 2] 内的所有根
  2. 使用多项式拟合指数函数 exe^x 在区间 [-1, 1] 上，比较不同次数的拟合效果
  3. 使用切比雪夫多项式逼近正弦函数
  4. 实现一个函数，使用多项式求函数在区间内的最大值
  5. 使用多项式积分计算 ∫01cos⁡(x)dx\int_0^1 \cos(x) dx 的近似值