python - 使用 SciPy 模块构建稀疏二维拉普拉斯矩阵

Question

我需要构建如下所示的 2D 拉普拉斯算子:

，哪里

, I 是单位矩阵。到目前为止，我已经使用 diags method of scipy 完成了它, 但我想知道是否有更聪明的方法来使用 block_diag method .有没有人试过用这种方法构建二维拉普拉斯算子？

我目前创建它的方法是通过这个函数:

from scipy.sparse import diags
# Defining the size of the matrix
nx = 3
ny = 3
N  = nx*ny
main_diag = [-4.0 for i in xrange(N)]
side_diag = []
for i in xrange(1,N):
    if i%4 == 0:
        side_diag.append(0)
    else:
        side_diag.append(1)
up_down_diag = [1.0 for i in xrange(N-4)]
diagonals = [main_diag,side_diag,side_diag,up_down_diag,up_down_diag]
laplacian = diags(diagonals, [0, -1, 1,nx,-nx], format="csr")
print laplacian.toarray()

Answer 1

我用数组替换了您对列表的使用:

import numpy as np
from scipy import sparse

nx, ny = 3, 3
N  = nx*ny
main_diag = np.ones(N)*-4.0
side_diag = np.ones(N-1)
side_diag[np.arange(1,N)%4==0] = 0
up_down_diag = np.ones(N-3)
diagonals = [main_diag,side_diag,side_diag,up_down_diag,up_down_diag]
laplacian = sparse.diags(diagonals, [0, -1, 1,nx,-nx], format="csr")
print laplacian.toarray()

生产

[[-4.  1.  0.  1.  0.  0.  0.  0.  0.]
 [ 1. -4.  1.  0.  1.  0.  0.  0.  0.]
 [ 0.  1. -4.  1.  0.  1.  0.  0.  0.]
 [ 1.  0.  1. -4.  0.  0.  1.  0.  0.]
 [ 0.  1.  0.  0. -4.  1.  0.  1.  0.]
 [ 0.  0.  1.  0.  1. -4.  1.  0.  1.]
 [ 0.  0.  0.  1.  0.  1. -4.  1.  0.]
 [ 0.  0.  0.  0.  1.  0.  1. -4.  0.]
 [ 0.  0.  0.  0.  0.  1.  0.  0. -4.]]

侧对角线的 [1 1 1 0 1 1 1 0] 模式是否正确？

对于像这样的小例子，它可能运行相同的速度，但使用数组而不是列表的大维度应该更快 - 而且它与 numpy 代码底层稀疏更一致。

像这样的统一对角线 diags 看起来很不错。

我只玩过另一个 SO 问题的 block 格式。 https://stackoverflow.com/a/34124377/901925

coo 适用于由重叠的较小矩阵组成的矩阵，例如有限元刚度。但是将对角线重铸为 coo 是乏味的。

就其值(value)而言，sparse.diags 使用 dia_matrix，将对角线列表转换为 dia data矩阵。你可以看看那个，但它的布局并不那么明显。为了制作一个 csr 矩阵，diags 将此 dia 格式转换为 coo，然后再转换为 csr 。但通常您不必担心所有这些转换。使用对您的问题最有意义的格式，并让 sparse 处理转换细节。

如果您想更多地探索 block 格式，您需要概述如何将您的问题视为 block 。