# Generalized linear algebraic equation solver


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. – Henrik Schumacher May 29 '18 at 22:56
• More detail on my problem here: math.stackexchange.com/questions/2801044/… – Kai May 29 '18 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 '18 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. – Daniel Lichtblau May 30 '18 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 '18 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[2]]
]


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 '18 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. – Henrik Schumacher May 31 '18 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 '18 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 '18 at 9:59