如何基于向量-矩阵映射实现函数计算与2×2子矩阵行列式求解

需要基于坐标计算函数f(x_i,y_i)的值，结果需呈现为如下矩阵形式：

f(x1,y1),f(x2,y1),...,f(xn,y1)
f(x1,y2),f(x2,y2),...,f(xn,y2)
..............................
f(x1,yn),f(x2,yn),...,f(xn,yn)

为避免循环计算，计划采用向量X=(x1,x2,x3,...,xn)和Y=(y1,y2,y3,...,yn)计算f(X,Y)，后续需基于每个f(xi,yi)的位置计算2×2子矩阵的行列式：

d_{ij}=f(x_i,y_j)+f(x_{i+1},y_{j+1})-f(x_i,y_{j+1})-f(x_{i+1},y_j)

并将结果存入新矩阵。梳理的步骤为：

• 获取向量X和Y；
• 计算f(X,Y)；
• 将f(X,Y)重新映射为矩阵结构；
• 计算每个行列式值。

已通过numpy.meshgrid完成第一步，但后续步骤实现存疑，给出的示例代码如下：

import numpy as np
import scipy as scp
def GG(r):
    return np.sqrt((y-y0)**2 + r**2 * np.sin(r))
def beta1(r):
    return ((y-y0) - r * np.cos(r)) / GG(r)
def f(x,y):
    return scp.integrate.quad_vec(lambda r: 1/(np.sinh(self.beta1(r))))
x = np.linspace(2,5,10)
y = np.linspace(-5,5,10)
y0 = 1
x0 = 1
a = np.stack(np.meshgrid(x, y, indexing="ij"), axis=-1)
xx = a[:, 0, 0]
yy = a[0, :, 1]
fxy = f(xx,yy)

第一步：修正网格矩阵生成 #

原代码提取xx和yy的方式错误，meshgrid配合indexing='ij'可直接生成对应所有(x_i,y_j)组合的二维矩阵，无需手动堆叠提取：

x = np.linspace(2,5,10)
y = np.linspace(-5,5,10)
y0 = 1
x0 = 1  # 原代码中x0未使用，若后续需要可保留

# 生成网格矩阵：X为(n,n)，每行都是x向量；Y为(n,n)，每列都是y向量
X, Y = np.meshgrid(x, y, indexing='ij')

第二步：修正函数定义（支持向量化计算） #

原代码函数存在三个核心问题：

1. GG和beta1依赖全局变量y，无法处理向量化输入；
2. f函数误用self（无类上下文），且quad_vec调用缺少积分上下限；
3. 未处理quad_vec的返回值结构（返回积分结果+误差）。

修正后的函数如下（假设积分区间为[0, 10]，可根据需求调整）：

import numpy as np
from scipy.integrate import quad_vec

def GG(r, y_val, y0):
    return np.sqrt((y_val - y0)**2 + r**2 * np.sin(r))

def beta1(r, y_val, y0):
    return ((y_val - y0) - r * np.cos(r)) / GG(r, y_val, y0)

def f(x_val, y_val, y0):
    # quad_vec返回(积分结果, 误差)，这里取积分结果
    result, _ = quad_vec(lambda r: 1/np.sinh(beta1(r, y_val, y0)), 0, 10)
    return result

第三步：计算完整的f(X,Y)矩阵 #

利用numpy向量化特性，直接将网格矩阵X和Y传入f函数，quad_vec支持向量输入，可一次性计算所有网格点的积分值：

# 计算所有(x_i,y_j)对应的f值，得到(n,n)的矩阵
f_matrix = f(X, Y, y0)

第四步：高效计算行列式d_ij矩阵 #

无需循环，利用numpy切片操作直接实现行列式公式，时间复杂度为O(n²)：

# 切片获取对应子矩阵的元素：
# f[:-1,:-1] 对应f(x_i,y_j)
# f[1:,1:]   对应f(x_{i+1},y_{j+1})
# f[:-1,1:]  对应f(x_i,y_{j+1})
# f[1:,:-1]  对应f(x_{i+1},y_j)
d_matrix = f_matrix[:-1, :-1] + f_matrix[1:, 1:] - f_matrix[:-1, 1:] - f_matrix[1:, :-1]

完整可运行代码 #

import numpy as np
from scipy.integrate import quad_vec

def GG(r, y_val, y0):
    return np.sqrt((y_val - y0)**2 + r**2 * np.sin(r))

def beta1(r, y_val, y0):
    return ((y_val - y0) - r * np.cos(r)) / GG(r, y_val, y0)

def f(x_val, y_val, y0):
    result, _ = quad_vec(lambda r: 1/np.sinh(beta1(r, y_val, y0)), 0, 10)
    return result

# 生成向量X和Y
x = np.linspace(2,5,10)
y = np.linspace(-5,5,10)
y0 = 1

# 生成网格矩阵
X, Y = np.meshgrid(x, y, indexing='ij')

# 计算f(X,Y)矩阵
f_matrix = f(X, Y, y0)

# 计算行列式d_ij矩阵
d_matrix = f_matrix[:-1, :-1] + f_matrix[1:, 1:] - f_matrix[:-1, 1:] - f_matrix[1:, :-1]

# 可选：打印结果验证
print("f矩阵形状:", f_matrix.shape)
print("d矩阵形状:", d_matrix.shape)

关键说明 #

• meshgrid(indexing='ij')确保生成的矩阵符合需求：X[i,j] = x_i，Y[i,j] = y_j，对应f(x_i,y_j)的位置；
• quad_vec支持向量化输入，可一次性计算所有网格点的积分，避免循环；
• 切片操作是numpy中高效处理矩阵局部元素的核心方式，完全替代循环计算行列式，大幅提升效率。

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