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

Question

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

，其中

，I为单位矩阵。到目前为止，我已经使用 diags method of scipy, but I wonder whether there is a smarter way to do it using the 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 问题的块格式。

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

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

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

Answer 2

N-dimensional拉普拉斯算子可以表示为一维拉普拉斯算子的克罗内克积：

import scipy.sparse as sp
def laplacian2D(N):
    diag=np.ones([N*N])
    mat=sp.spdiags([diag,-2*diag,diag],[-1,0,1],N,N)
    I=sp.eye(N)
    return sp.kron(I,mat,format='csr')+sp.kron(mat,I)

Connectivity pattern of 2D Laplacian - img

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

Buildin a sparse 2D laplacian matrix using SciPy modules

python

numpy

scipy

sparse-matrix