使用对角线覆盖构建 Scipy 矩阵

我想在偏移位置覆盖几个相同 4x4 大小的数组，创建一个大数组，其中重叠的元素添加在一起。

[a] 0 0 0
 0 [a+b] 0 0
 0 0 [b+c] 0
 0 0 0 [c+...

这是为了创建用于连续梁分析的刚度矩阵。使用循环创建这很容易，但是使用 Scipy 稀疏矩阵例程是否有更好的方法？ block_diag 很接近，但没有提供覆盖的方法。

Members = [[100.0, 1000.0],[100.0, 2000.0],[200.0, 2000.0]]

SSM = np.zeros((2*M+2,2*M+2), dtype = float)
i=0
for mem in Members:
    a = SMX(mem)
    print a
    print "*"*40

    SSM[i:i+4,i:i+4] = SSM[i:i+4,i:i+4] + a
    print SSM
    i +=2

def SMX(member):
    """ Prismatic beam member stiffness function
    Purpose: Create a member stiffness matrix for a given set of member parameters.
    Input:  Modulus of elasticity - E
            Moment of inertia with respect to the z axis - Iz
            Length of the span - L

    Returns: 4 x 4 matrix representing the member stiffness with respect to end conditions    
"""

    L = member[0]
    Iz = member[1]

    pbms = np.zeros((4,4), dtype = float)

    pbms[0,0] = (12.0 * E * Iz)/L**2
    pbms[0,1] = (6.0 * E * Iz)/L**2
    pbms[0,2] = -pbms[0,0]
    pbms[0,3] = pbms[0,1]

    pbms[1,0] = pbms[0,1]
    pbms[1,1] = (4.0 * E * Iz)/L
    pbms[1,2] = -pbms[0,1]
    pbms[1,3] = (2.0 * E * Iz)/L

    pbms[2,0] = pbms[0,2]
    pbms[2,1] = -pbms[0,1]
    pbms[2,2] = pbms[0,0]
    pbms[2,3] = -pbms[0,1]

    pbms[3,0] = pbms[0,1]
    pbms[3,1] = pbms[1,3]
    pbms[3,2] = -pbms[0,1]
    pbms[3,3] = pbms[1,1]

    return pbms

结果输出：

[[  12000.    6000.  -12000.    6000.]
 [   6000.  400000.   -6000.  200000.]
 [ -12000.   -6000.   12000.   -6000.]
 [   6000.  200000.   -6000.  400000.]]
****************************************
[[  12000.    6000.  -12000.    6000.       0.       0.       0.       0.]
 [   6000.  400000.   -6000.  200000.       0.       0.       0.       0.]
 [ -12000.   -6000.   12000.   -6000.       0.       0.       0.       0.]
 [   6000.  200000.   -6000.  400000.       0.       0.       0.       0.]
 [      0.       0.       0.       0.       0.       0.       0.       0.]
 [      0.       0.       0.       0.       0.       0.       0.       0.]
 [      0.       0.       0.       0.       0.       0.       0.       0.]
 [      0.       0.       0.       0.       0.       0.       0.       0.]]
[[  24000.   12000.  -24000.   12000.]
 [  12000.  800000.  -12000.  400000.]
 [ -24000.  -12000.   24000.  -12000.]
 [  12000.  400000.  -12000.  800000.]]
****************************************
[[   12000.     6000.   -12000.     6000.        0.        0.        0.
         0.]
 [    6000.   400000.    -6000.   200000.        0.        0.        0.
         0.]
 [  -12000.    -6000.    36000.     6000.   -24000.    12000.        0.
         0.]
 [    6000.   200000.     6000.  1200000.   -12000.   400000.        0.
         0.]
 [       0.        0.   -24000.   -12000.    24000.   -12000.        0.
         0.]
 [       0.        0.    12000.   400000.   -12000.   800000.        0.
         0.]
 [       0.        0.        0.        0.        0.        0.        0.
         0.]
 [       0.        0.        0.        0.        0.        0.        0.
         0.]]
[[   6000.    3000.   -6000.    3000.]
 [   3000.  400000.   -3000.  200000.]
 [  -6000.   -3000.    6000.   -3000.]
 [   3000.  200000.   -3000.  400000.]]
****************************************
[[   12000.     6000.   -12000.     6000.        0.        0.        0.
         0.]
 [    6000.   400000.    -6000.   200000.        0.        0.        0.
         0.]
 [  -12000.    -6000.    36000.     6000.   -24000.    12000.        0.
         0.]
 [    6000.   200000.     6000.  1200000.   -12000.   400000.        0.
         0.]
 [       0.        0.   -24000.   -12000.    30000.    -9000.    -6000.
      3000.]
 [       0.        0.    12000.   400000.    -9000.  1200000.    -3000.
    200000.]
 [       0.        0.        0.        0.    -6000.    -3000.     6000.
     -3000.]
 [       0.        0.        0.        0.     3000.   200000.    -3000.
    400000.]]

刚刚注意到我遗漏了上面的 M = 3 。

我没有看到问题scipy.sparse.block_diag http://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.block_diag.html，但我可能不太明白你所说的“覆盖”是什么意思。例如，您在问题中说明的输出可以这样形成：

In [41]: from scipy import sparse

In [42]: a = np.array([[1, 2], [3, 4]])

In [43]: b = np.array([[5, 6], [7, 8]])

In [44]: c = np.array([[9, 10], [11, 12]])

In [45]: m = sparse.block_diag([a, a+b, b+c])

In [46]: m.A
Out[46]: 
array([[ 1,  2,  0,  0,  0,  0],
       [ 3,  4,  0,  0,  0,  0],
       [ 0,  0,  6,  8,  0,  0],
       [ 0,  0, 10, 12,  0,  0],
       [ 0,  0,  0,  0, 14, 16],
       [ 0,  0,  0,  0, 18, 20]])

Edit：浏览了一些有关光束分析的材料后，我想我更好地理解了您想要如何覆盖矩阵。

例如，假设有四个 2x2 数组 a、b、c 和 d。 “覆盖”的示例如下：

a[0,0]     a[0,1]            0              0            0
a[1,0]  a[1,1]+b[0,0]     b[0,1]            0            0
  0        b[1,0]      b[1,1]+c[0,0]     c[0,1]          0
  0           0           c[1,0]      c[1,1]+d[0,0]  d[0,1]
  0           0              0           d[1,0]      d[1,1]

形成这个的一种方法是block_diag在这里进行了演示。它涉及一个表达式，每个 2x2 数组都有一个术语，可以将其重写为循环，因此它并不能真正回答您的问题。

In [92]: a,b,c,d
Out[92]: 
(array([[1, 2],
       [3, 4]]),
 array([[5, 6],
       [7, 8]]),
 array([[ 9, 10],
       [11, 12]]),
 array([[13, 14],
       [15, 16]]))

In [93]: bd = sparse.block_diag  # A convenient alias

In [94]: m = bd([a,0,0,0])+bd([0,b,0,0])+bd([0,0,c,0])+bd([0,0,0,d])

In [95]: m.A
Out[95]: 
array([[ 1,  2,  0,  0,  0],
       [ 3,  9,  6,  0,  0],
       [ 0,  7, 17, 10,  0],
       [ 0,  0, 11, 25, 14],
       [ 0,  0,  0, 15, 16]])