Block tridiagonal matrices, how to program this kind of matrix?

Question

We consider h=1/(n+1) and n>=1. We consider a matrix K such as :

and the matrix L = 1/h² Id, where Id is the identity matrix.

We would like to program the matrix A :

I can program K and L, but how to program A? A is a matrix composed of matrices. How to implement this in Python?

Possible duplicate of [Block tridiagonal matrix python](https://stackoverflow.com/questions/5842903/block-tridiagonal-matrix-python) — RagingRoosevelt, Dec 11 '17 at 20:58
@RagingRoosevelt It doesn't look like a dupe to me. This is asking for a genuine block matrix while linked OP doesn't seem to know what a block matrix actually is. — Paul Panzer, Dec 11 '17 at 21:33
Yeah, looks like you're right. I think I may have misread OP. — RagingRoosevelt, Dec 11 '17 at 21:35
scipy.sparse now has a `block_diag` function: https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.block_diag.html. You may need to create a `K -L/-L K` subblock first. sparse block assembly uses the `coo` sparse format to define a new sparse matrix. — hpaulj, Dec 12 '17 at 00:46

Paul Panzer · Accepted Answer · 2017-12-11T22:07:48.523

If you are open to numpy/scipy (recommended):

Use scipy.linalg.toeplitz to create banded matrices and numpy.kron to create patterns of repeated blocks:

import numpy as np
import scipy.linalg

h = 10
K = np.zeros((4,))
K[:2] = (4 / h**2, -1 / h**2)

K = scipy.linalg.toeplitz(K)

L = np.identity(4) / h**2

KK = np.identity(3)

LL = scipy.linalg.toeplitz((0, -1, 0))

A = np.kron(LL, L) + np.kron(KK, K)

# array([[ 0.04, -0.01,  0.  ,  0.  , -0.01, -0.  ,  0.  ,  0.  ,  0.  ,  0.  ,  0.  ,  0.  ],
#        [-0.01,  0.04, -0.01,  0.  , -0.  , -0.01, -0.  ,  0.  ,  0.  ,  0.  ,  0.  ,  0.  ],
#        [ 0.  , -0.01,  0.04, -0.01,  0.  , -0.  , -0.01, -0.  ,  0.  ,  0.  ,  0.  ,  0.  ],
#        [ 0.  ,  0.  , -0.01,  0.04,  0.  ,  0.  , -0.  , -0.01,  0.  ,  0.  ,  0.  ,  0.  ],
#        [-0.01, -0.  ,  0.  ,  0.  ,  0.04, -0.01,  0.  ,  0.  , -0.01, -0.  ,  0.  ,  0.  ],
#        [-0.  , -0.01, -0.  ,  0.  , -0.01,  0.04, -0.01,  0.  , -0.  , -0.01, -0.  ,  0.  ],
#        [ 0.  , -0.  , -0.01, -0.  ,  0.  , -0.01,  0.04, -0.01,  0.  , -0.  , -0.01, -0.  ],
#        [ 0.  ,  0.  , -0.  , -0.01,  0.  ,  0.  , -0.01,  0.04,  0.  ,  0.  , -0.  , -0.01],
#        [ 0.  ,  0.  ,  0.  ,  0.  , -0.01, -0.  ,  0.  ,  0.  ,  0.04, -0.01,  0.  ,  0.  ],
#        [ 0.  ,  0.  ,  0.  ,  0.  , -0.  , -0.01, -0.  ,  0.  , -0.01,  0.04, -0.01,  0.  ],
#        [ 0.  ,  0.  ,  0.  ,  0.  ,  0.  , -0.  , -0.01, -0.  ,  0.  , -0.01,  0.04, -0.01],
#        [ 0.  ,  0.  ,  0.  ,  0.  ,  0.  ,  0.  , -0.  , -0.01,  0.  ,  0.  , -0.01,  0.04]])

If it has to be pure Python:

Make a matrix of matrices and the use zip to transpose the middle dimensions and flat list comprehension to make 2D.

K = [[0.04 if i==j else -0.01 if i-j in {-1, 1} else 0.0 for j in range(4)] for i in range(4)]

L = [[-0.01 if i==j else 0.0 for j in range(4)] for i in range(4)]

Z = [[0.0 for j in range(4)] for i in range(4)]

# matrix of matrices
A = [[K if i==j else L if i-j in {-1, 1} else Z for j in range(3)] for i in range(3)]

# make 2d
A = [[a_IJij for a_IJi in a_Ii for a_IJij in a_IJi] for a_I in A for a_Ii in zip(*a_I)]

for a in A:
    print(' '.join(f'{i:5.2f}' for i in a))

#  0.04 -0.01  0.00  0.00 -0.01  0.00  0.00  0.00  0.00  0.00  0.00  0.00
# -0.01  0.04 -0.01  0.00  0.00 -0.01  0.00  0.00  0.00  0.00  0.00  0.00
#  0.00 -0.01  0.04 -0.01  0.00  0.00 -0.01  0.00  0.00  0.00  0.00  0.00
#  0.00  0.00 -0.01  0.04  0.00  0.00  0.00 -0.01  0.00  0.00  0.00  0.00
# -0.01  0.00  0.00  0.00  0.04 -0.01  0.00  0.00 -0.01  0.00  0.00  0.00
#  0.00 -0.01  0.00  0.00 -0.01  0.04 -0.01  0.00  0.00 -0.01  0.00  0.00
#  0.00  0.00 -0.01  0.00  0.00 -0.01  0.04 -0.01  0.00  0.00 -0.01  0.00
#  0.00  0.00  0.00 -0.01  0.00  0.00 -0.01  0.04  0.00  0.00  0.00 -0.01
#  0.00  0.00  0.00  0.00 -0.01  0.00  0.00  0.00  0.04 -0.01  0.00  0.00
#  0.00  0.00  0.00  0.00  0.00 -0.01  0.00  0.00 -0.01  0.04 -0.01  0.00
#  0.00  0.00  0.00  0.00  0.00  0.00 -0.01  0.00  0.00 -0.01  0.04 -0.01
#  0.00  0.00  0.00  0.00  0.00  0.00  0.00 -0.01  0.00  0.00 -0.01  0.04

score 1 · Answer 2 · answered Dec 11 '17 at 21:04

Use the numpy.eye function to get an identity matrix with ones on the diagonal and zeros everywhere else. You can use the optional parameter k=x to specify the offset for the diagonal. If you multiply these identity matrices by a scalar, you'll be able to introduce the L and K. Add the matrices together to get the tri-diagonal matrix that you're looking for.

import numpy as np

m, n = 4, 4
L = 1
K = 4
A = np.eye(m, n, k=-1) * (-L) + np.eye(m, n) * K + np.eye(m, n, k=1) * (-L)

Block tridiagonal matrices, how to program this kind of matrix?

2 Answers2