如何拆分Numpy刚度矩阵以实现静力凝聚？

嘿，这事儿用Numpy的索引切片就能轻松搞定，完全不用手动敲子矩阵的元素，既高效又不会出错！我来给你一步步演示：

首先，你的原始刚度矩阵是5x5的，根据你给出的拆分需求，其实就是把矩阵按「前3行/列（保留自由度rr）」和「后2行/列（凝聚自由度cc）」来划分。直接用Numpy的切片语法就能快速提取出各个子矩阵：

完整实现代码 #

from sympy import *
import numpy as np

# 你原来的定义代码
u1 = Symbol("u1")
u2 = Symbol("u2")
q3 = Symbol("q3")
u4 = Symbol("u4")
q5 = Symbol("q5")

disp_vector = np.array([u1, u2, q3, u4, q5])

stiffness_matrix = np.array([[1, 0, 0, -1, 0],
                             [0, 0.12, 0.6, 0, 0.6],
                             [0, 0.6, 4, 0, 2],
                             [-1, 0, 0, 1, 0],
                             [0, 0.6, 2, 0, 4]])

force_vector = np.array([0, 40, -26.7, 0, 0])

# 核心：用切片拆分刚度矩阵
krr = stiffness_matrix[:3, :3]  # 前3行，前3列
krc = stiffness_matrix[:3, 3:]  # 前3行，后2列
kcr = stiffness_matrix[3:, :3]  # 后2行，前3列
kcc = stiffness_matrix[3:, 3:]  # 后2行，后2列

# 可以打印出来验证，和你手动定义的完全一致！
print("krr:\n", krr)
print("\nkrc:\n", krc)
print("\nkcr:\n", kcr)
print("\nkcc:\n", kcc)

灵活扩展：如果自由度不是连续的怎么办？ #

要是你的保留自由度和凝聚自由度不是像这样前后连续的（比如分散的索引，比如r=[0,2,4]，c=[1,3]），可以用np.ix_()来生成索引网格，这样能正确提取子矩阵：

# 举个例子，假设保留自由度是索引0、2、4，凝聚是1、3
r = [0, 2, 4]
c = [1, 3]
krr = stiffness_matrix[np.ix_(r, r)]
krc = stiffness_matrix[np.ix_(r, c)]
kcr = stiffness_matrix[np.ix_(c, r)]
kcc = stiffness_matrix[np.ix_(c, c)]

这样不管自由度怎么分布，都能准确拆分出需要的子矩阵，完美适配静力凝聚的需求~

备注：内容来源于stack exchange，提问作者Hür Doğan ÜNLÜ