有限元法门式框架分析：矩阵逆与向量相乘维度不匹配问题

我正在开发采用有限元法（FEM）的门式框架分析代码，在尝试将矩阵特定位置的逆与数组相乘时，Python返回错误：shapes (9,9) and (9,1,1) not aligned: 9 (dim 1) != 1 (dim 1)。

我有一个15×15的矩阵kportico和一个15行的数组f，但仅需对索引[3:11]位置的元素进行相乘操作。已确认矩阵和数组内容正确，但无法定位问题所在。相关代码片段如下：

k = [k1global, k2global, k3global, k4global]

kportico = np.zeros([15, 15])

for i, m in enumerate(k):
    kportico[i*3:i*3+6, i*3:i*3+6] += m

f = np.array([[f1global[0]],
              [f1global[1]],
              [f1global[2]],
              [f1global[3]+f2globalx[0]+f2globaly[0]],
              [f1global[4]+f2globalx[1]+f2globaly[1]],
              [f1global[5]+f2globalx[2]+f2globaly[2]],
              [f2globalx[3]+f2globaly[3]+f3globalx[0]+f3globaly[0]],
              [f2globalx[4]+f2globaly[4]+f3globalx[1]+f3globaly[1]],
              [f2globalx[5]+f2globaly[5]+f3globalx[2]+f2globaly[2]],
              [f4global[3]+f3globalx[3]+f3globaly[3]],
              [f4global[4]+f3globalx[4]+f3globaly[4]],
              [f4global[5]+f3globalx[5]+f3globaly[5]],
              [f4global[0]],
              [f4global[1]],
              [f4global[2]]])

u = np.dot(np.linalg.pinv(kportico[np.ix_([3, 4, 5, 6, 7, 8, 9, 10, 11], [3, 4, 5, 6, 7, 8, 9, 10, 11])]), f[np.ix_([3, 4, 5, 6, 7, 8, 9, 10, 11])])

问题原因 #

报错核心是数组维度不匹配：

• 切片后的子矩阵是(9,9)的二维数组
• 用np.ix_处理f后得到的是(9,1,1)的三维数组，矩阵乘法要求第二个参数的第一维度等于第一个参数的第二维度（即9），但这里是1，因此触发维度对齐错误。

修复方法 #

有两种简单的方式解决维度问题：

方法1：直接切片代替np.ix_ #

对于列向量f，直接用f[3:12]（索引3到11共9个元素）即可得到(9,1)的二维数组，与子矩阵维度完美匹配：

# 提取9×9的子矩阵
k_sub = kportico[3:12, 3:12]
# 提取9×1的子数组
f_sub = f[3:12]
# 计算位移u
u = np.dot(np.linalg.pinv(k_sub), f_sub)

方法2：压缩多余维度 #

如果一定要保留np.ix_的使用，可以用np.squeeze()去掉多余的维度，将(9,1,1)转为(9,1)：

# 提取子矩阵
k_sub = kportico[np.ix_([3,4,5,6,7,8,9,10,11], [3,4,5,6,7,8,9,10,11])]
# 提取子数组并压缩维度
f_sub = np.squeeze(f[np.ix_([3,4,5,6,7,8,9,10,11])])
# 或者用reshape强制调整维度
# f_sub = f[np.ix_([3,4,5,6,7,8,9,10,11])].reshape(9,1)
# 计算位移u
u = np.dot(np.linalg.pinv(k_sub), f_sub)

补充说明 #

np.ix_的设计初衷是生成多维索引以实现矩阵的不规则切片，当仅对一维列向量使用时，会额外增加一个维度，导致数组从二维变为三维，完全没必要。直接切片f[3:12]更简洁高效，完全满足需求。

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