Трехдиагональная Matrix с матричными блоками

Question

Я прочитал this, но я не был в состоянии создать (N^2 х N^2) - матрица с (п х N) - матрицЯ на нижняя и верхняя боковые диагональные и T по диагонали. Я пробовал это

def prep_matrix(N): 
    I_N = np.identity(N) 
    NV = zeros(N*N - 1) 
    # now I want NV[0]=NV[1]=...=NV[N-1]:=I_N 

но я не знаю, как заполнить NV моими матрицами. Что я могу сделать? Я много узнал о том, как создавать тридиагональные матрицы со скалярами, но не с матричными блоками.

Answer 1

Вы можете инициализировать так:

n = 3 
I, T, A = np.identity(n), np.ones(n), np.zeros([n * n, n * n]) 
for i in range(n): 
    a, b, c = i * n, (i + 1) * n, (i + 2) * n 
    A[a:b, a:b] = T 
    if i < n - 1: A[b:c, a:b] = A[a:b, b:c] = I 

В приведенном выше примере имеет следующий вывод:

[[ 1. 1. 1. 1. 0. 0. 0. 0. 0.] 
[ 1. 1. 1. 0. 1. 0. 0. 0. 0.] 
[ 1. 1. 1. 0. 0. 1. 0. 0. 0.] 
[ 1. 0. 0. 1. 1. 1. 1. 0. 0.] 
[ 0. 1. 0. 1. 1. 1. 0. 1. 0.] 
[ 0. 0. 1. 1. 1. 1. 0. 0. 1.] 
[ 0. 0. 0. 1. 0. 0. 1. 1. 1.] 
[ 0. 0. 0. 0. 1. 0. 1. 1. 1.] 
[ 0. 0. 0. 0. 0. 1. 1. 1. 1.]] 

Answer 2

Ради искусства, другой подход, основанный на bmat():

Z = np.zeros((n, n)) 
I = np.eye(n) 
T = np.ones((n, n))*2 

B = np.diag([1]*n) + np.diag([2]*np.ones(n-1), 1) + np.diag(2*np.ones(n-1), -1) 
# B = array([[ 1., 2., 0.], 
#   [ 2., 1., 2.], 
#   [ 0., 2., 1.]]) 
# build 2d list `Bm` replacing 0->Z, 1->T, 2->I: 
bdict = {0.:Z, 1.:T, 2.:I} 
Bm = [[bdict[i] for i in rw] for rw in B] 
# Use the power of bmat to construct matrix: 
A = np.asarray(np.bmat(Bm))