块三对角矩阵，如何编程这种矩阵？答案

【问题标题】：Block tridiagonal matrices, how to program this kind of matrix?块三对角矩阵，如何编程这种矩阵？
【发布时间】：2018-05-25 11:00:28
【问题描述】：

我们考虑h=1/(n+1) 和n>=1。我们考虑一个矩阵K，例如：

还有矩阵L = 1/h² Id，其中Id是单位矩阵。

我们想编程矩阵A：

我可以编程K 和L，但是如何编程A？ A是由矩阵组成的矩阵。如何在 Python 中实现？

【问题讨论】：

Block tridiagonal matrix python的可能重复
@RagingRoosevelt 这对我来说不像是骗子。这是要求一个真正的块矩阵，而链接的 OP 似乎并不知道块矩阵实际上是什么。
是的，看起来你是对的。我想我可能误读了 OP。
scipy.sparse 现在有一个block_diag 函数：docs.scipy.org/doc/scipy/reference/generated/…。您可能需要先创建一个K -L/-L K 子块。稀疏块组装使用coo稀疏格式来定义一个新的稀疏矩阵。

标签： python numpy math matrix

【解决方案1】：

如果您对numpy/scipy（推荐）持开放态度：

使用scipy.linalg.toeplitz 创建带状矩阵，使用numpy.kron 创建重复块的模式：

import numpy as np
import scipy.linalg

h = 10
K = np.zeros((4,))
K[:2] = (4 / h**2, -1 / h**2)

K = scipy.linalg.toeplitz(K)

L = np.identity(4) / h**2

KK = np.identity(3)

LL = scipy.linalg.toeplitz((0, -1, 0))

A = np.kron(LL, L) + np.kron(KK, K)

# array([[ 0.04, -0.01,  0.  ,  0.  , -0.01, -0.  ,  0.  ,  0.  ,  0.  ,  0.  ,  0.  ,  0.  ],
#        [-0.01,  0.04, -0.01,  0.  , -0.  , -0.01, -0.  ,  0.  ,  0.  ,  0.  ,  0.  ,  0.  ],
#        [ 0.  , -0.01,  0.04, -0.01,  0.  , -0.  , -0.01, -0.  ,  0.  ,  0.  ,  0.  ,  0.  ],
#        [ 0.  ,  0.  , -0.01,  0.04,  0.  ,  0.  , -0.  , -0.01,  0.  ,  0.  ,  0.  ,  0.  ],
#        [-0.01, -0.  ,  0.  ,  0.  ,  0.04, -0.01,  0.  ,  0.  , -0.01, -0.  ,  0.  ,  0.  ],
#        [-0.  , -0.01, -0.  ,  0.  , -0.01,  0.04, -0.01,  0.  , -0.  , -0.01, -0.  ,  0.  ],
#        [ 0.  , -0.  , -0.01, -0.  ,  0.  , -0.01,  0.04, -0.01,  0.  , -0.  , -0.01, -0.  ],
#        [ 0.  ,  0.  , -0.  , -0.01,  0.  ,  0.  , -0.01,  0.04,  0.  ,  0.  , -0.  , -0.01],
#        [ 0.  ,  0.  ,  0.  ,  0.  , -0.01, -0.  ,  0.  ,  0.  ,  0.04, -0.01,  0.  ,  0.  ],
#        [ 0.  ,  0.  ,  0.  ,  0.  , -0.  , -0.01, -0.  ,  0.  , -0.01,  0.04, -0.01,  0.  ],
#        [ 0.  ,  0.  ,  0.  ,  0.  ,  0.  , -0.  , -0.01, -0.  ,  0.  , -0.01,  0.04, -0.01],
#        [ 0.  ,  0.  ,  0.  ,  0.  ,  0.  ,  0.  , -0.  , -0.01,  0.  ,  0.  , -0.01,  0.04]])

如果它必须是纯 Python：

制作矩阵矩阵并使用zip转置中间维度和平面列表理解以制作二维。

K = [[0.04 if i==j else -0.01 if i-j in {-1, 1} else 0.0 for j in range(4)] for i in range(4)]

L = [[-0.01 if i==j else 0.0 for j in range(4)] for i in range(4)]

Z = [[0.0 for j in range(4)] for i in range(4)]

# matrix of matrices
A = [[K if i==j else L if i-j in {-1, 1} else Z for j in range(3)] for i in range(3)]

# make 2d
A = [[a_IJij for a_IJi in a_Ii for a_IJij in a_IJi] for a_I in A for a_Ii in zip(*a_I)]

for a in A:
    print(' '.join(f'{i:5.2f}' for i in a))

#  0.04 -0.01  0.00  0.00 -0.01  0.00  0.00  0.00  0.00  0.00  0.00  0.00
# -0.01  0.04 -0.01  0.00  0.00 -0.01  0.00  0.00  0.00  0.00  0.00  0.00
#  0.00 -0.01  0.04 -0.01  0.00  0.00 -0.01  0.00  0.00  0.00  0.00  0.00
#  0.00  0.00 -0.01  0.04  0.00  0.00  0.00 -0.01  0.00  0.00  0.00  0.00
# -0.01  0.00  0.00  0.00  0.04 -0.01  0.00  0.00 -0.01  0.00  0.00  0.00
#  0.00 -0.01  0.00  0.00 -0.01  0.04 -0.01  0.00  0.00 -0.01  0.00  0.00
#  0.00  0.00 -0.01  0.00  0.00 -0.01  0.04 -0.01  0.00  0.00 -0.01  0.00
#  0.00  0.00  0.00 -0.01  0.00  0.00 -0.01  0.04  0.00  0.00  0.00 -0.01
#  0.00  0.00  0.00  0.00 -0.01  0.00  0.00  0.00  0.04 -0.01  0.00  0.00
#  0.00  0.00  0.00  0.00  0.00 -0.01  0.00  0.00 -0.01  0.04 -0.01  0.00
#  0.00  0.00  0.00  0.00  0.00  0.00 -0.01  0.00  0.00 -0.01  0.04 -0.01
#  0.00  0.00  0.00  0.00  0.00  0.00  0.00 -0.01  0.00  0.00 -0.01  0.04

【讨论】：

【解决方案2】：

使用numpy.eye 函数获得一个单位矩阵，其中对角线为 1，其他位置为 0。您可以使用可选参数k=x 指定对角线的偏移量。如果将这些单位矩阵乘以一个标量，您将能够引入L 和K。将这些矩阵相加即可得到您要查找的三对角矩阵。

import numpy as np

m, n = 4, 4
L = 1
K = 4
A = np.eye(m, n, k=-1) * (-L) + np.eye(m, n) * K + np.eye(m, n, k=1) * (-L)

【讨论】：