如何使用循环在矩阵中编码矩阵？

Question

所以我使用两个变量 u（向量，x - y）和标量 k 创建了这个 3x3 G 矩阵（此处未显示，它与我的问题无关）。 x_j = (x_1 (j), x_2 (j), x_3 (j)) 和 y_j = (y_1 (j), y_2 (j), y_3 (j))。 alpha_j 是一个 3x3 矩阵。 A 矩阵是大小为 3nx3n 的块对角矩阵。我在使用 W 矩阵时遇到问题。如何编写大小为 3nx3n 的矩阵，其中第 (i,j) 个块是由 alpha_i*G_[ij]*alpha_j 给出的 3x3 矩阵？我迷路了。

我的 alpha_j 矩阵似乎也有问题。循环不断向我抛出错误，“只有长度为 1 的数组可以转换为 Python 标量。”请帮助：/

def W(x, y, k, alpha, A):
    u = x - y
    n = x.shape[0]
    W = np.zeros((3*n, 3*n))

    for i in range(0, n-1):
        for j in range(0, n-1):
                #u = -np.array([[x[i,0] - x[j,0]], [x[i,1] -               x[j,1]], [0]]) ??
                W[i][j] = (alpha_j(alpha, A) * G(u, k) * alpha_j(alpha, A))
        W[i][i] = np.zeros((n, n))

    return W

def alpha_j(a, A):
    alph = np.array([[0,0,0],[0,0,0],[0,0,0]],complex)
    rho = np.random.rand(3,1)
    for i in range(0, 2):
        for j in range(0, 2):
            alph[i][j] = (rho[i] * a * A[i][j])
    return alph

#-------------------------------------------------------------------

x1 = np.array([[1], [2], [0]])
y1 = np.array([[4], [5], [0]])

# SYSTEM PARAMETERS

# incoming Wave angle
theta = 0 # can range from [0, 2pi)

# susceptibility
chi = 10 + 1j

# wavelength
lam = 0.5 # microns (values between .4-.7)

# frequency
k = (2 * np.pi)/lam # 1/microns

# volume
V_0 = (0.05)**3 # microns^3

# incoming wave vector
K = k * np.array([[0], [np.sin(theta)], [np.cos(theta)]])

# polarization vector
vecinc = np.array([[1], [0], [0]]) # (can choose any vector perpendicular to K)

# for the fixed alpha case
alpha = (V_0 * 3 * chi)/(chi + 3)

# 3 x 3 matrix
A = np.matlib.identity(3) # could be any symmetric matrix, 


#-------------------------------------------------------------------

# TEST FUNCTIONS

test = G((x1-y1), k)
print(test)

w = W(x1, y1, k, alpha, A)
print(w)

有时我的 W 循环会向我抛出错误，“无法设置具有序列的数组元素”。但是我需要将这个任意矩阵 W 中的每个数组元素设置为通过将 alpha 乘以 G 创建的 3x3 矩阵...

Answer 1

对于如何为每个元素创建一个带有块的新数组的问题，以下应该可以解决问题：

G = np.random.random([3,3])
result = np.zeros([9,9])
num_blocks = 3
a = np.random.random([3,3])
b = np.random.random([3,3])
for i in range(G.shape[0]):
    for j in range(G.shape[1]):
        block_result = a*G[i,j]*b
        for k in range(num_blocks):
            for l in range(num_blocks):
                result[3*i + k, 3*j + l] = block_result[i, j]

您应该能够从那里进行概括。我希望我理解正确。

编辑：看来我没有正确理解。我留下它，希望它会刺激你的答案。总体思路是生成要操作的索引范围，然后直接对它们进行操作。切片也可能有帮助。

啊，你问如何创建一个填充块的对角线。在这种情况下：

num_diagonal_blocks = 3 # for example

for block_dim in range(num_diagonal_blocks)
    # do your block calculation...
    for k in range(G.shape[0]):
        for l in range(G.shape[1]):
            result[3*block_dim + k, 3*block_dim + l] = # assign to element of block

我想差不多了。