Legendre多项式牛顿法求根报错：无法广播数组形状(11,10)到(11)

看起来你的代码有两个核心问题：一是Legendre多项式的递推逻辑不完整，二是返回了不必要的多维数组导致后续迭代中x的形状混乱，最终触发广播错误。我来一步步帮你修正：

错误根源分析 #

1. Legendre函数的递推逻辑缺失：你当前只初始化了P0、P1，然后直接跳步计算PN，但中间的P2到PN-1都是初始的0值，这会导致PN的计算完全错误。
2. 不必要的多维数组返回：你返回了存储所有阶数多项式的二维数组，但牛顿法只需要N阶的函数值和导数值。这会让后续的negfn = -a/b得到二维数组，和一维的x相加后，x变成二维数组，第二次调用Legendre时就会触发形状不匹配的错误。

修正后的完整代码 #

1. 修正Legendre多项式函数 #

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

def Legendre(N, x):
    # 初始化0阶和1阶Legendre多项式及其导数
    P_prev_prev = 1 + 0*x  # P₀(x) = 1
    P_prev = x             # P₁(x) = x
    Q_prev_prev = 0*x      # P₀的导数为0
    Q_prev = 1 + 0*x       # P₁的导数为1
    
    # 处理低阶特殊情况
    if N == 0:
        return P_prev_prev, Q_prev_prev
    if N == 1:
        return P_prev, Q_prev
    
    # 递推计算从2阶到N阶的多项式和导数
    for k in range(2, N+1):
        P_current = ((2*k - 1)*x*P_prev - (k - 1)*P_prev_prev) / k
        Q_current = (2*k - 1)*P_prev + Q_prev_prev
        # 更新前序项，准备下一轮迭代
        P_prev_prev, P_prev = P_prev, P_current
        Q_prev_prev, Q_prev = Q_prev, Q_current
    
    # 仅返回N阶的函数值和导数值（与输入x形状一致）
    return P_prev, Q_prev

2. 修正牛顿法求解函数 #

def quadrature(N, tol=1e-19):
    # 初始猜测：使用切比雪夫根作为Legendre根的近似，收敛更快
    i = np.arange(N)
    x = np.cos(np.pi*(4*i + 3)/(4*N + 2))
    
    eps = np.inf
    iteration = 0
    imax = 20
    
    while eps > tol and iteration < imax:
        # 获取N阶Legendre函数值和导数值
        P_N, Q_N = Legendre(N, x)
        # 牛顿迭代公式：x = x - f(x)/f'(x)
        delta = -P_N / Q_N
        x_new = x + delta
        # 计算迭代的最大误差，判断是否收敛
        eps = np.max(np.abs(delta))
        x = x_new
        iteration += 1
        # 可选：打印迭代过程
        # print(f"Iteration {iteration}, max error: {eps:.2e}")
    
    # 计算高斯-勒让德积分权重
    _, Q_N = Legendre(N, x)
    w = 2 / ((1 - x**2) * Q_N**2)
    return x, w

关键修正点说明 #

• Legendre函数：
  • 去掉了存储所有阶数的二维数组，只保留当前和前两个阶数的结果，既节省内存又保证递推逻辑正确。
  • 新增循环完整计算从2阶到N阶的多项式和导数，确保结果准确。
  • 仅返回与输入x形状一致的一维数组，避免后续广播错误。
• quadrature函数：
  • 明确了牛顿迭代的数学公式，确保迭代逻辑清晰。
  • 新增了误差计算逻辑，能正确判断迭代是否收敛（原代码中eps未更新，会导致循环逻辑失效）。
  • 权重计算部分也适配了修正后的Legendre函数返回值。

现在调用quadrature(10)就能正常返回10阶高斯-勒让德积分的根和权重了。

内容的提问来源于stack exchange，提问作者john sam