Мне нужно построить конкретную матрицу, где некоторые элементы сами являются матрицами [python].Блочная трехдиагональная матрица

Question

У меня есть та часть кода, которая генерирует B и I:

import numpy as np
import scipy.sparse.linalg as sc

size = 10

#1. create the matrix B and I
def matrix(size):
    dx = size / (size - 1)
    B = np.identity(size ) * (-4)
    counter = 0
    sizeCounter = size
    for i in range(1, size):
        B[i, counter] = 1
        B[counter, i] = 1
        counter += 1
    B[0,:]=0
    B[:,0]=0
    B[size-1,:]=0
    B[:,size-1]=0

    I=np.identity(size)
    I[0,:]=0
    I[:,0]=0
    I[size-1,:]=0
    I[:,size-1]=0


    return B,I

Я хотел бы сгенерировать матрицу M в квадрате (размер ^ 2) x (размер ^ 2). Так что, если размер равен 10М будет 100х100.М будет иметь следующую форму на картинке. .Как я могу это сделать?

Answer 1

Вот три варианта:

Сначала создайте строительные блоки:

Я использую einsum, чтобы написать диагонали.

>>> import numpy as np
>>> 
>>> size = 6
>>> 
>>> O, B, I = OBI = np.zeros((3, size, size))
>>> np.einsum('ii->i', I[1:-1, 1:-1])[:] = 1
>>> np.einsum('ii->i', B[1:-1, 1:-1])[:] = -4
>>> np.einsum('ii->i', B[1:-2, 2:-1])[:] = 1
>>> np.einsum('ii->i', B[2:-1, 1:-2])[:] = 1

Далее,построить контейнер:

Я использую toeplitz для построения (скалярной) трехдиагональной матрицы.

>>> from scipy import linalg
>>>
>>> T = np.zeros((size,), dtype=int)
>>> T[:2] = 1, 2
>>> T = linalg.toeplitz(T)

Решение 1:

Объединить с помощью np.block:

>>> M1 = np.block(list(map(list, OBI[T])))

Решение 2:

Объединение с использованием reshape:

>>> M2 = OBI[T].swapaxes(1, 2).reshape(size*size, size*size)

Решение 3:

Использование np.kron:

>>> Tm1 = np.zeros((size,))
>>> Tm1[1] = 1
>>> Tm1 = linalg.toeplitz(Tm1)
>>> 
>>> M3 = np.kron(np.identity(size,), B) + np.kron(Tm1, I)

Answer 2

Вы можете использовать np.concatenate для построения каждой строки вашей матрицы матриц.

B, I = matrix(size)
A = np.zeros((size, size))

M = []
for i in range(size):
    if i == 0:
        M.append(np.concatenate([ B, I, np.tile(A, (1, size-2)) ], axis = -1 ))
    elif i == size-1:
        M.append(np.concatenate([ np.tile(A, (1, size-2)), I, B ], axis = -1))
    else:
        M.append(np.concatenate([ np.tile(A, (1, i-1)), I, B, I, np.tile(A, (1, size-i-2))], axis = -1))
M = np.concatenate(M, axis = 0)

Использование списка:

M = np.concatenate([
        np.concatenate([ B, I, np.tile(A, (1, size-2)) ], axis = -1 ) if i == 0 else
        np.concatenate([ np.tile(A, (1, size-2)), I, B ], axis = -1) if i == size-1
        else np.concatenate([ np.tile(A, (1, i-1)), I, B, I, np.tile(A, (1, size-i-2))], axis = -1)
        for i in range(size)], axis = 0)