如何用NumPy无循环实现二维Gauss-Legendre积分点与权重计算

现有带循环的实现代码 #

import numpy as np

def get_quad_points_weights(n):
    # 获取一维Gauss-Legendre点和权重
    x, w = np.polynomial.legendre.leggauss(n)
    quad_points = np.zeros((n*n, 2))
    quad_weights = np.zeros(n*n)
    idx = 0
    for i in range(n):
        for j in range(n):
            quad_points[idx, 0] = x[i]
            quad_points[idx, 1] = x[j]
            quad_weights[idx] = w[i] * w[j]
            idx += 1
    return quad_points, quad_weights

# n=2时的输出示例
points_2, weights_2 = get_quad_points_weights(2)
print("n=2时的积分点：")
print(points_2)
print("n=2时的权重：")
print(weights_2)

n=2时的预期输出 #

n=2时的积分点：
[[-0.57735027 -0.57735027]
[-0.57735027 0.57735027]
[ 0.57735027 -0.57735027]
[ 0.57735027 0.57735027]]
n=2时的权重：
[1. 1. 1. 1.]

常见错误的无循环尝试（供参考） #

def get_quad_points_weights_no_loop_wrong(n):
    x, w = np.polynomial.legendre.leggauss(n)
    quad_points = np.array([x, x]).reshape(-1, 2)  # 错误：未生成所有点的组合
    quad_weights = np.array([w * w]).flatten()     # 错误：未计算所有权重的乘积组合
    return quad_points, quad_weights

正确的无循环实现方法 #

利用NumPy的网格生成和外积操作，完全替代循环逻辑，代码如下：

import numpy as np

def get_quad_points_weights_no_loop(n):
    x, w = np.polynomial.legendre.leggauss(n)
    # 生成与原循环顺序一致的二维网格点（indexing='ij'保证行优先）
    x_grid, y_grid = np.meshgrid(x, x, indexing='ij')
    # 将网格点展平并组合为(n², 2)的积分点数组
    quad_points = np.column_stack((x_grid.flatten(), y_grid.flatten()))
    # 计算权重的外积并展平，得到所有w[i]*w[j]的组合
    quad_weights = np.outer(w, w).flatten()
    return quad_points, quad_weights

# 测试n=2的情况
points_2_no_loop, weights_2_no_loop = get_quad_points_weights_no_loop(2)
print("无循环实现n=2时的积分点：")
print(points_2_no_loop)
print("无循环实现n=2时的权重：")
print(weights_2_no_loop)

关键逻辑说明 #

• meshgrid(x, x, indexing='ij')：生成和原双层循环顺序完全匹配的网格点，避免因索引顺序导致的点排列混乱。
• np.column_stack：将两个一维网格数组拼接为二维积分点数组，格式与循环实现一致。
• np.outer(w, w)：直接生成权重的乘积矩阵，展平后即为所有权重组合，等价于循环中w[i]*w[j]的计算。

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