我想像这样生成分区对角矩阵A
给出矩阵B
B = -np.diag(np.ones(n - 2), -1) - np.diag(np.ones(n - 2), 1) + 4 * np.diag(np.ones(n - 1))
例如,
有没有一种方法可以不使用循环?
抱歉,第一次错误地上传了矩阵A和B的图形.
解决方法:
您可以将构建块堆叠到查找表中,然后通过在其中建立索引来构建A:
>>> from scipy import sparse
>>>
>>> n = 5
>>> B = sparse.diags([-1, 4, -1], [-1, 0, 1], (n-1, n-1), dtype=int).A
>>> A = sparse.diags([1, 2, 1], [-1, 0, 1], (n-1, n-1), dtype=int).A
# 0 means 0 0 0 ...,
# 1 means -I
# 2 means B
>>>
# next line builds the lookup table (using np.stack)
# does the lookup ...[A]
# and flattens the resulting 4D array after swapping
# the middle axes; the swap reorders the entries from
# Vert, Horz, vert, horz to Vert, vert, Horz, horz
>>> A = np.stack([np.zeros_like(B), -np.identity(n-1, int), B])[A].swapaxes(1, 2).reshape((n-1)*(n-1), -1)
>>> A
array([[ 4, -1, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[-1, 4, -1, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[ 0, -1, 4, -1, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, 0],
[ 0, 0, -1, 4, 0, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0, 0],
[-1, 0, 0, 0, 4, -1, 0, 0, -1, 0, 0, 0, 0, 0, 0, 0],
[ 0, -1, 0, 0, -1, 4, -1, 0, 0, -1, 0, 0, 0, 0, 0, 0],
[ 0, 0, -1, 0, 0, -1, 4, -1, 0, 0, -1, 0, 0, 0, 0, 0],
[ 0, 0, 0, -1, 0, 0, -1, 4, 0, 0, 0, -1, 0, 0, 0, 0],
[ 0, 0, 0, 0, -1, 0, 0, 0, 4, -1, 0, 0, -1, 0, 0, 0],
[ 0, 0, 0, 0, 0, -1, 0, 0, -1, 4, -1, 0, 0, -1, 0, 0],
[ 0, 0, 0, 0, 0, 0, -1, 0, 0, -1, 4, -1, 0, 0, -1, 0],
[ 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, -1, 4, 0, 0, 0, -1],
[ 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, 0, 4, -1, 0, 0],
[ 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, -1, 4, -1, 0],
[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, -1, 4, -1],
[ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1, 0, 0, -1, 4]])
请注意,稀疏构造函数仅出于其便利性而使用.稀疏矩阵将立即转换为密集矩阵(使用.A属性).
标签:python,numpy