使用Scipy.sparse.diags实现Scipy三对角矩阵

Question

我在使用numpy数组创建三对角矩阵时遇到了困难。我已经成功地复制了这里提供的结果，但是我无法将这些技术应用到我的问题上。我可能也误解了scipy.sparse.diag的应用。

为了背景，我正在解决一个问题，需要生成三对角矩阵，以使用有限差分数值方法来数值求解常微分方程。

from scipy.sparse import diags
import numpy as np

v1 = [3*i**2 +(i/2) for i in range(1, 6)]
v2 = [-(6*i**2 - 1) for i in range(1, 6)]
v3 = [3*i**2 -(i/2) for i in range(1, 6)]

matrix = np.array([v1, v2, v3])

matrix 等于什么。

array([[3.5,   13. ,   28.5,   50. ,   77.5],
       [-5. ,  -23. ,  -53. ,  -95. , -149. ],
       [2.5,   11. ,   25.5,   46. ,   72.5]])

浏览了Scipy文档和上述链接中的示例后，我原以为以下代码会产生Tridiagonal_1，但实际得到的是Tridiagonal_2。

diags(matrix, [-1,0,1], (5, 5)).toarray() 

期望的Tridiagonal_1：

array([[  -5. ,    2.5 ,     0. ,    0. ,     0. ],
       [  13. ,   -23. ,    11. ,    0. ,     0. ],
       [   0. ,    28.5.,  -53. ,   25.5,     0. ],
       [   0. ,    0. ,     50 ,   -95.,     46. ],
       [   0. ,    0. ,      0. ,   77.5., -149. ]])

代码产生了Tridiagonal_2：

array([[  -5. ,    2.5,    0. ,    0. ,    0. ],
       [   3.5,  -23. ,   11. ,    0. ,    0. ],
       [   0. ,   13. ,  -53. ,   25.5,    0. ],
       [   0. ,    0. ,   28.5,  -95. ,   46. ],
       [   0. ,    0. ,    0. ,   50. , -149. ]])

我原本期望 offset = [-1,0,1] 将对角线上的元素向左移动，但第一个偏移量将第一个 diag 移动到了下一行。这是正确的吗？还是我的代码有误导致了这种行为？

- Josmoor98

---

diags将输入视为数组列表，并在不移位的情况下使用它们。而spdiags则将2D数组映射到对角线上，顶部会有溢出。这只是处理长度不同的对角线的不同约定。 - hpaulj

Answer 1

你所得到的输出似乎符合文档，特别是文档中的示例所说的。 你可以使用spdiags来获得你想要的结果：

from scipy import sparse

matrix = np.array([[3.5,   13. ,   28.5,   50. ,   77.5],
                   [-5. ,  -23. ,  -53. ,  -95. , -149. ],
                   [2.5,   11. ,   25.5,   46. ,   72.5]]

sparse.spdiags(matrix, (1,0,-1), 5, 5).T.A
# array([[  -5. ,    2.5,    0. ,    0. ,    0. ],
#        [  13. ,  -23. ,   11. ,    0. ,    0. ],
#        [   0. ,   28.5,  -53. ,   25.5,    0. ],
#        [   0. ,    0. ,   50. ,  -95. ,   46. ],
#        [   0. ,    0. ,    0. ,   77.5, -149. ]])

- Paul Panzer