# Generalized linear algebraic equation solver

$\newcommand{\d}{\vec{d}}$ $\newcommand{\S}{\vec{S}}$

In Mathematica one can easily solve a linear system given by $$A \vec{S} = \vec{d}$$ where $A$ is a matrix, simply by using S = LinearSolve[A,d].

I don't know if this has a name, but I am wondering if there is a nice way to solve a generalized system of linear equations such as this one (this is a specific example but the form of the matrix on the left is unimportant, only that it is a "matrix of matrices"): $$\begin{bmatrix} 1 & A & 0 & & & 0 \\[8pt] -A & 1 & A & & & \\[1pt] 0 & -A & 1 & \ddots & & \\[1pt] & & \ddots & \ddots & & \\[1pt] & & & & 1 & A \\[6pt] 0 & & & & -A & 1 \end{bmatrix} \begin{bmatrix} \S_1 \\[3pt] \S_2 \\[3pt] \S_3 \\[3pt] \vdots \\[3pt] \S_{N-1} \\[3pt] \S_{N} \end{bmatrix} = \begin{bmatrix} \d_1 \\[3pt] \d_2 \\[3pt] \d_3 \\[3pt] \vdots \\[3pt] \d_{N-1} \\[3pt] \d_{N} \end{bmatrix}$$ where $1$ is the identity matrix and $A$ is an $n\times n$ matrix. The $\S_i$ and $\d_i$ are $n$-component vectors, the $\d_i$ are known and I want to solve for the $\S_i$.

The "brute force" way to solve this is to Flatten the "vectors of vectors" and ArrayFlatten the "matrix of matrices", then use LinearSolve, then use Partition[S,n] to recover the vectors. But I feel like such systems must not be that uncommon and maybe there is some built in method of doing this.

• By the way, this does not look like Crank-Nicolson. Crank-Nicolson is for transient PDEs but your system looks more like stemming from an elliptic boundary value problem. May 29, 2018 at 22:56
• More detail on my problem here: math.stackexchange.com/questions/2801044/…
– Kai
May 29, 2018 at 23:34
• I am simply following an assignment which says to solve using Crank-Nicolson, I believe that my code is working, the results do seem reasonable. I applied boundary conditions that $\rho$ and $v$ are zero at the edges.
– Kai
May 29, 2018 at 23:35
• Some actual code would be useful here. Otherwise it becomes a do-my-homework type of question, better suited for math.SE perhaps. May 30, 2018 at 15:34
• I already solved the problem using a method somewhat similar to the one below. My question is asking if such a generalized linear algebra problem can be solved in Mathematica without using Flatten and Partition. It would be nice to have an option for such a problem built in to LinearSolve for example.
– Kai
May 30, 2018 at 22:38

You can build the system matrix as follows

A = RandomReal[{-1, 1}, {2, 2}];
A = A\[Transpose].A;
n = 10;
AA = Plus[
KroneckerProduct[SparseArray[{Band[{2, 1}] -> -1, Band[{1, 2}] -> 1}, {n, n}], A],
KroneckerProduct[IdentityMatrix[n, SparseArray], N@IdentityMatrix]
]


The right-hand side can be assembled as follows:

d = RandomReal[{-1, 1}, {n, 2}];
S1 = RandomReal[{-1, 1}, {2}];
Sn = RandomReal[{-1, 1}, {2}];
b = Flatten[d];
b[[;; 2]] += A.S1;
b[[-2 ;;]] -= A.Sn;


Solve as usual (this is a usual linear system). I'd suggest a solver specialized for banded matrices:

x = LinearSolve[AA, b, Method -> "Banded"];


In order to convert back to a list of pairs, you can use

Partition[x,2]

• Thanks, while this does solve the problem (I wrote a similar code) I was hoping there might be a built in way to handle such problems.
– Kai
May 30, 2018 at 22:41
• More general block matrices can be converted to a a single large matrix with ArrayFlatten. In the end, such a task has to be done since LinearSolve uses external libraries which expect a matrix as input. May 31, 2018 at 9:54
• Aha I have not really used ArrayFlatten, but that makes the problem relatively trivial to do then. Just Flatten/ArrayFlatten everything, use LinearSolve, then Partition. Okay that's not too bad, especially since the matrix doesn't need to be partitioned again. I wonder if such a problem has a name, I can't imagine this is uncommon.
– Kai
May 31, 2018 at 9:57
• Sorry for the hassle but I changed the question to make it more general, so parts of your answer don't exactly make sense now in the context of the question, if you have the time to update it.
– Kai
May 31, 2018 at 9:59