求助：使用线性优化构造N次多项式逼近cos(x)的Python实现

嘿，我看你正在尝试用线性规划来构造N次多项式逼近cos(x)，这个思路挺靠谱的！先帮你把问题理得更清晰，然后把代码补全并修正可能的问题~

首先回顾你的优化问题：我们要找一个N次多项式 $P(x) = a_0 + a_1x + ... + a_Nx^N$，使得在[0,1]区间的M个采样点$x_i$上，$P(x_i)$和$\cos(x_i)$的误差绝对值被$e_i$覆盖，目标是最小化所有$e_i$的和。对应的线性规划形式是：

$$
\min \sum_{i=1}^{M} e_i
$$
约束条件：

1. $P(x_i) + e_i \geq \cos(x_i)$ → 等价于 $-P(x_i) - e_i \leq -\cos(x_i)$
2. $P(x_i) - e_i \leq \cos(x_i)$

接下来咱们把这个问题转换成scipy.optimize.linprog能处理的标准形式（linprog要求目标是最小化线性函数，约束为$A_{ub} \cdot x \leq b_{ub}$）：

变量定义 #

我们的优化变量是一个长度为$(N+1)+M$的数组：

• 前$N+1$个元素：多项式系数$[a_0, a_1, ..., a_N]$
• 后$M$个元素：误差变量$[e_1, e_2, ..., e_M]$

目标函数系数 #

目标是最小化$\sum e_i$，所以目标系数数组c的前$N+1$位是0，后$M$位是1。

约束矩阵与右侧向量 #

我们需要为每个采样点生成2个约束，总共$2M$个约束：

• 对于第$i$个采样点的第一个约束（$-P(x_i) - e_i \leq -\cos(x_i)$）：对应矩阵的第$i$行，多项式系数位置是$-x_i^0, -x_i^1, ..., -x_i^N$，第$N+1+i$位（对应$e_i$）是-1，其余位置是0
• 对于第$i$个采样点的第二个约束（$P(x_i) - e_i \leq \cos(x_i)$）：对应矩阵的第$M+i$行，多项式系数位置是$x_i^0, x_i^1, ..., x_i^N$，第$N+1+i$位是-1，其余位置是0

完整Python实现代码 #

from scipy.optimize import linprog
import numpy as np
import matplotlib.pyplot as plt

def polynomial_approx_cos(M, N, interval=[0, 1]):
    # 生成M个均匀采样点
    x_points = np.linspace(interval[0], interval[1], M)
    cos_vals = np.cos(x_points)
    
    # 定义变量数量：N+1个系数 + M个误差项
    num_vars = N + 1 + M
    
    # 目标函数系数：最小化e_i的和
    c = np.zeros(num_vars)
    c[N+1:] = 1  # 后M个变量是e_i，系数为1
    
    # 构造约束矩阵A_ub和右侧向量b_ub
    A_ub = np.zeros((2*M, num_vars))
    b_ub = np.zeros(2*M)
    
    for i in range(M):
        x = x_points[i]
        # 第一个约束：-a0 -a1x - ... -aNx^N - ei <= -cos(xi)
        A_ub[i, :N+1] = -np.power(x, np.arange(N+1))
        A_ub[i, N+1 + i] = -1
        b_ub[i] = -cos_vals[i]
        
        # 第二个约束：a0 +a1x + ... +aNx^N - ei <= cos(xi)
        A_ub[M+i, :N+1] = np.power(x, np.arange(N+1))
        A_ub[M+i, N+1 + i] = -1
        b_ub[M+i] = cos_vals[i]
    
    # 所有变量非负？其实e_i必须非负，系数a可以是任意实数，所以只约束e_i >=0
    # 构造边界：a的边界是(-inf, inf)，e的边界是(0, inf)
    bounds = [(-np.inf, np.inf)]*(N+1) + [(0, np.inf)]*M
    
    # 求解线性规划
    result = linprog(c, A_ub=A_ub, b_ub=b_ub, bounds=bounds, method='highs')
    
    if result.success:
        # 提取多项式系数
        coeffs = result.x[:N+1]
        # 构造多项式函数（polyval是从高次到低次，所以翻转系数）
        def P(x):
            return np.polyval(np.flip(coeffs), x)
        
        # 可视化结果
        x_plot = np.linspace(interval[0], interval[1], 100)
        plt.plot(x_plot, np.cos(x_plot), label='cos(x)')
        plt.plot(x_plot, P(x_plot), label=f'{N}-degree polynomial')
        plt.scatter(x_points, cos_vals, color='red', marker='x', label='Sample points')
        plt.legend()
        plt.xlabel('x')
        plt.ylabel('y')
        plt.title(f'Linear Programming Approximation of cos(x) with {N}-degree polynomial')
        plt.show()
        
        return coeffs, result.x[N+1:]
    else:
        print("Optimization failed:", result.message)
        return None, None

# 测试：用M=20个采样点，构造3次多项式逼近
coeffs, errors = polynomial_approx_cos(M=20, N=3)
if coeffs is not None:
    print("Polynomial coefficients (a0, a1, a2, a3):", coeffs)
    print("Maximum error:", np.max(errors))

代码说明 #

1. 采样点生成：用linspace在[0,1]生成均匀分布的M个点，也可以换成随机点
2. 约束构造：循环每个采样点，分别生成两个约束的行向量和右侧值
3. 变量边界：误差项e_i必须非负（因为是绝对值误差的上界），多项式系数a可以是任意实数，所以设置对应的边界
4. 求解器选择：使用method='highs'，这是scipy较新版本推荐的高效线性规划求解器
5. 可视化：把原函数、逼近多项式和采样点画在一起，直观查看逼近效果

你可以调整M（采样点数量）和N（多项式次数）来观察不同的逼近效果，比如增大N会让多项式更贴合cos(x)，增大M会让逼近在更多点上满足误差约束~

备注：内容来源于stack exchange，提问作者Felipe Oliveira