This auxiliary section reminds the reader when a system of linear algebraic equations has a solution. Such systems are called consistent.

R solves linear systems

A system of m linear equations in n unknowns \( x_1 , \ x_2 , \ \ldots , \ x_n \) is a set of m equations of the form

\begin{equation} \label{EqSolve.1} \begin{cases} a_{11} \,x_1 + a_{12}\, x_2 + \cdots + a_{1n}\, x_n &= b_1 , \\ a_{21} \,x_1 + a_{22}\, x_2 + \cdots + a_{2n}\, x_n &= b_2 , \\ \qquad \vdots & \qquad \vdots \\ a_{m1} \,x_1 + a_{m2}\, x_2 + \cdots + a_{mn}\, x_n &= b_m , \end{cases} \end{equation}

The numbers \( a_{ij} \) are known as the coefficients of the system. The m×n matrix \( {\bf A} = [\,a_{ij}\,] , \) whose \( (i,\, j) \) entry is the coefficient \( a_{ij} \) of the system of linear equations is called the coefficient matrix, and is denoted by

\[ {\bf A} = \left[ \begin{array}{cccc} a_{1,1} & a_{1,2} & \cdots & a_{1,n} \\ a_{2,1} & a_{2,2} & \cdots & a_{2,n} \\ \vdots & \vdots & \ddots & \vdots \\ a_{m,1} & a_{m,2} & \cdots & a_{m,n} \end{array} \right] \qquad\mbox{or} \qquad {\bf A} = \left( \begin{array}{cccc} a_{1,1} & a_{1,2} & \cdots & a_{1,n} \\ a_{2,1} & a_{2,2} & \cdots & a_{2,n} \\ \vdots & \vdots & \ddots & \vdots \\ a_{m,1} & a_{m,2} & \cdots & a_{m,n} \end{array} \right) . \]

Let \( {\bf x} =\left[ x_1 , x_2 , \ldots x_n \right]^{\mathrm T} \) be the column vector of unknowns, where "T" stands for transposition. Then the product \( {\bf A}\,{\bf x} \) of the \( m \times n \) coefficient matrix A and the \( n \times 1 \) column vector x is an \( m \times 1 \) matrix (column-vector)

\[ \left[ \begin{array}{cccc} a_{1,1} & a_{1,2} & \cdots & a_{1,n} \\ a_{2,1} & a_{2,2} & \cdots & a_{2,n} \\ \vdots & \vdots & \ddots & \vdots \\ a_{m,1} & a_{m,2} & \cdots & a_{m,n} \end{array} \right] \, \begin{bmatrix} x_1 \\ x_2 \\ \vdots \\ x_n \end{bmatrix} = \begin{bmatrix} a_{1,1} x_1 + a_{1,2} x_2 + \cdots + a_{1,n} x_n \\ a_{2,1} x_1 + a_{2,2} x_2 + \cdots + a_{2,n} x_n \\ \vdots \\ a_{m,1} x_1 + a_{m,2} x_2 + \cdots + a_{m,n} x_n \end{bmatrix} , \]

whose entries are the right-hand sides of our system of linear equations.

If we define another column vector \( {\bf b} = \left[ b_1 , b_2 , \ldots b_m \right]^{\mathrm T} \) whose components are the right-hand sides \( b_{i} ,\) the given system is equivalent to the vector equation

\begin{equation} \label{EqSolve.2} {\bf A} \, {\bf x} = {\bf b} . \end{equation}

We say that \( s_1 , \ s_2 , \ \ldots , \ s_n \) is a solution of the system if all \( m \) equations hold true when

\[ x_1 = s_1 , \ x_2 = s_2 , \ \ldots , \ x_n = s_n . \]

Sometimes a system of linear equations is known as a set of simultaneous equations; such terminology emphasizes that a solution is an assignment of values to each of the \( n \) unknowns such that each and every equation holds with this assignment. It is also referred to simply as a linear system.

Mathematica has a built-in command LinearSolve[A, b] that solves a linear system of equations, given a coefficient matrix A and a vector of constants b. We need to learn some more theory before we can entirely understand this command, but we can begin to explore its use. For now, consider these methods experimental and do not let it replace row-reducing augmented matrices:

\[ {\bf A}\, {\bf x} = {\bf b} , \]

where A is a given matrix and b is specified column vector. The command LinearSolve[A, b] will provide the vector x of solutions to the linear system \( {\bf A}\,{\bf x} = {\bf b} \) of equations with coefficient matrix A and the vector of constants b.

The system of algebraic equations Ax = b may have a unique solution x, may have infinite many solutions, or may have no solution at all.

Example 1: Let us take a look at the system of algebraic equations

\[ \begin{cases} \phantom{2\,}x+2\,y &=3 , \\ 2\,x + 4\,y &=6 . \end{cases} \tag{1.1} \]

The two equations define the same line on the x-y plane, meaning that we have infinitely many solutions x = 3−2y. The corresponding singular (its determinant is zero) matrix A and nonhomogeneous term b are

\[ {\bf A} = \begin{bmatrix} 1 & 2 \\ 2 & 4 \end{bmatrix} , \qquad {\bf b} = \begin{bmatrix} 3 \\ 6 \end{bmatrix} . \]

Let's apply the Solve command to this system of algebraic equations. However, the command Solve is within a special pacage "matlib"

install.packages("matlib", dependencies = TRUE)
A <- matrix(c(1, 2, 2, 4), 2, 2)
b <- c(3, 6)
Solve(A, b, fractions = TRUE)

Now we apply the Solve command to a system with no solutions:

x1 = LinearSolve[A, b]

Out[8]= {3, 0}

So we see that the vector v₁ = [ 3, 0 ] gives a solution to system (1.1) of equations. Is this solution unique? To answer this question, we need to find the kernel (or null space) of matrix A:

NullSpace[A]

Out[10]= {{-2, 1}}

Hence, the null space is spanned on the vector [ -2, 1 ]^T and the general solution to the given system of linear equations Ax = b is

\[ x = 3\,t , \quad y = t , \]

where t is an arbitrary real number. We can also find the general solution with the following commands:

A = {{1, 2}, {2, 4}}
b = {3, 6}
Solve[A.{x, y} == b]

However, if we take another vector b = [−3, 6]^T. Obviously, the new system of equations

\[ \begin{cases} \phantom{2\,}x+2\,y &=-3 , \\ 2\,x + 4\,y &=\phantom{-}6 . \end{cases} \tag{1.2} \]

defines two parallel lines on x-y-plane that do not intersect, meaning there is no solut. Therefore, system (1.2) has no solution. We check whether Mathematica has a capability to determine this conclusion.

b2 = {-3, 6}
x2 = LinearSolve[A, b2]

Out[10]=
LinearSolve::nosol: Linear equation encountered that has no solution.

Therefore, the system has no solution for this vector b₂ because matrix rank

MatrixRank[A]

Out[12]= 1

is 1, which does not match the dimensions 2×2 of the given matrix. ■

Example 2: Let us consider the singular matrix \( {\bf A} = \begin{bmatrix} 2&\phantom{-}1&2 \\ 0&\phantom{-}1&0 \\ 1&-2& 1 \end{bmatrix} .\) We want to determine when the linear equation \( {\bf A} \, {\bf x} = {\bf b} \) has a solution; namely, we need to find conditions on the vector b so the equation \( {\bf A} \, {\bf x} = {\bf b} \) has a solution. We use R to answer this question.

The command rank provides the rank of a matrix A. We will need the latter command when matrix A would have complex entries.

x = NullSpace[Atranspose]

Out[6]= {{-1, 5, 2}}

So we get the null space of the transpose matrix to be spanned on the vector [-1, 5, 2 ]. The system of linear equations Ax = b has a solution (in this case, we say the the given system is consistent) if and only if the vector b is orthogonal to every vector from the kernel of the transpose (in general adjoint) matrix. In other words, its dot product with each such vector must be zero. So we ask Mathematica to do this job:

b = {b1, b2, b3};
Dot[x, y] == 0

Out[8]= {-b1 + 5 b2 + 2 b3} == 0

x.y == 0

Of course, we can achieve the same result by solving the system of linear equations Ax = b directly using Gaussian elimination method. Mathematica does not provide this algorithmitically fastest way to solve a linear algebraic equation; instead it uses Gauss--Jordan elimination procedure, which is more computationally demanded (and practically is not used). In this case, Mathematica offers a special command to achieve this: RowReduce[A]. ■

A linear system of algebraic equations is said to be consistent if it has at least one solution and inconsistent if it has no solution.

Thus, a consistent linear system has either one solution or infinitely many solutions---there are no other possibilities.

Theorem 1: Let A be an m × n matrix.

(Overdetermined case) If m > n, then the linear system \( {\bf A}\,{\bf x} = {\bf b} \) is inconsistent for at least one vector b in \( \mathbb{R}^n . \)
(Undetermined case) If m < n, then for each vector b in \( \mathbb{R}^m \) the linear system \( {\bf A}\,{\bf x} = {\bf b} \) is either inconsistent or has infinite many solutions.

Theorem 2: If A is an \( m \times n \) matrix with \( m < n, \) then \( {\bf A}\,{\bf x} = {\bf 0} \) has infinitely many solutions.

Theorem 3: A linear system of algebraic equations \( {\bf A}\,{\bf x} = {\bf b} \) is consistent if and only if the vector b is orthogonal to every solution y of the adjoint homogeneous equation A*y = 0, where A* is the adjoint matrix to A. If A is real, then the adjoint matrix is just transposed and the statement is equivalent to \( {\bf A}^{\mathrm T} \,{\bf y} = {\bf 0} \) or \( {\bf y}^T {\bf A} = {\bf 0} . \)

Theorem 3 provides us an algorithm of determination whether a system of equations is consistent or not.

The kernel of a matrix A (which considered as a linear map), also known as the null spacenullspace, is the set of vectors v such that Av = 0 (so matrix A maps into zero).
The cokernel of a matrix A is the left null space, of a matrix A that consists of all column vectors y such that A*y = 0 or y^TA = 0^T when A is real, where T denotes the transpose of a matrix. The left null space of A is the same as the kernel of A* (or A^T if A is real).

Example 3: We consider a singular 3×3 matrix and 3-column vector b

\[ {\bf A} = \begin{bmatrix} 1&\phantom{-}1&1 \\ 2&\phantom{-}3&1 \\ 0&-1&1 \end{bmatrix} \qquad \mbox{and} \qquad {\bf b} = \begin{bmatrix} 3\\ 4 \\ 2 \end{bmatrix} . \]

We calculated its rank (number of linearly independent either rows or columns) So we see that the kernel of A^T is spanned on the vector

\[ {\bf y} = \left[ -2, 1, 1 \right]^{\mathrm T} . \]

However, if we take another input vector b = [3, 4, 2]^T, the corresponding system becomes consistent. ■

Example 4: Let us consider the matrix A and vector b:

\[ {\bf A} = \begin{bmatrix} \phantom{-}1&-2&3&-4 \\ -1&\phantom{-}1&1&\phantom{-}2 \\ \phantom{-}3&-5&5&-10 \\ \phantom{-}1 &-4&11&-8 \end{bmatrix} \qquad\mbox{amd} \qquad {\bf b} = \begin{bmatrix} 1 \\ 2 \\ 3 \\ 4 \end{bmatrix} . \]

The matrix A has rank 2.

So we see that the kernel of A^T is spanned on the vectors

\[ {\bf y}_1 = \left[ -3, -2, 0, 1 \right]^{\mathrm T} \qquad\mbox{and}\qquad {\bf y}_2 = \left[ -2, 1, 1, 0 \right]^{\mathrm T} . \]

End of Example 4

■

Theorem 4: The general solution x of a consistent linear system of algebraic equations \( {\bf A}\,{\bf x} = {\bf b} \) can be expressed as the sum of a particular solution x₀ of that system and the general solution of the corresponding homogeneous system

\[ {\bf x} = {\bf x}_0 + c_1 {\bf v}_1 + c_2 {\bf v}_2 + \cdots + c_k {\bf v}_k . \]

Here vector \( \{ {\bf v}_1 , \ldots , {\bf v}_k \} \) is a basis for the kernel (= null space) of A and \( c_1 , \ldots , c_k \) are arbitrary constants.

Theorem 5: Let A be an \( m \times n \) matrix. A linear system of equations \( {\bf A}\,{\bf x} = {\bf b} \) is consistent if and only if the rank of the coefficient matrix A is the same as the rank of the augmented matrix \( \left[ {\bf A}\,\vert \, {\bf b} \right] . \) ▣

The actual use of the term augmented matrix appears to have been introduced by the American mathematician Maxime Bôcher (1867--1918) in 1907.

Example 5: The linear system

\[ \begin{cases} x_1 + x_2 + x_3 + x_4 + x_5 &= 3 , \\ 2\,x_1 + x_2 + x_3 + x_4 + 2\, x_5 &= 4 , \\ x_1 - x_2 - x_3 + x_4 + x_5 &= 5 , \\ x_1 + x_4 + x_5 &= 4 , \end{cases} \tag{3.1} \]

is an example of a system of four equations in five unknowns, \( x_1 , \ x_2 , \ x_3 , \ x_4 , \ x_5 . \) One solution of this system is

\[ x_1 = -1 , \ x_2 = -2 , \ x_3 =1, \ x_4 = 3 , \ x_5 = 2 , \]

as you can easily verify by substituting these values into the equations. Every equation is satisfied for these values of \( x_1 , \ x_2 , \ x_3 , \ x_4 , \ x_5 . \) However, there is not the only solution to this system of equations. There are many others.

On the other hand, the system of linear equations

\[ \begin{cases} x_1 + x_2 + x_3 + x_4 + x_5 &= 3 , \\ 2\,x_1 + x_2 + x_3 + x_4 + 2\, x_5 &= 4 , \\ x_1 - x_2 - x_3 + x_4 + x_5 &= 5 , \\ x_1 + x_4 + x_5 &= 6 , \end{cases} \tag{3.2} \]

has no solution. There are no numbers we can assign to the unknowns \( x_1 , \ x_2 , \ x_3 , \ x_4 , \ x_5 \) so that all four equations are satisfied. We can solve these two equations simultaneously by considering the augmented matrix. An efficient way to do this is to form the partitioned matrix \( \left[ {\bf A} \, \vert \, {\bf b}_1 \, \vert {\bf b}_2 \right] \) in which the coefficient matrix A is "augmented" by all 2 vectors. Then this augmented matrix is reduced using Gauss--Jordan elimination to

\[ \left( \begin{array}{ccccc|cc} 1&1&1&1&1&3&3 \\ 2&1&1&1&2&4&4 \\ 1&-1&-1&1&1&5&5 \\ 1&0&0&1&1&4&6 \end{array} \right) \quad \Longrightarrow \quad \left( \begin{array}{ccccc|cc} 1&0&0&0&1&1&0 \\ 0&1&1&0&0&-1&0 \\ 0&0&0&1&0&3&0 \\ 0&0&0&0&0&0&1 \end{array} \right) . \]

The sixth column gives the required solution for the former system of equations. The last row in reduced form shows that the latter system of equations has no solution. ■

Equations in blocks

Suppose that we have the linear system \( {\bf M}\,{\bf z} = {\bf r} , \) in which the given coefficient matrix M is partitioned into 2 by 2 blocks: \( {\bf M} = \begin{bmatrix} {\bf A} & {\bf B} \\ {\bf C} & {\bf D} \end{bmatrix} , \) where A and D are square invertible matrices. Suppose that the right-hand side given vector is also broken as r = [ p q ]^T, while column vector z = [ x y ]^T is unknown. The system is equivalent to the pair of vector equations:

\begin{equation} \label{EqSolve.3} \begin{split} {\bf A}\, {\bf x} + {\bf B}\, {\bf y} = {\bf p} , \\ {\bf C}\, {\bf x} + {\bf D}\, {\bf y} = {\bf q} . \end{split} \end{equation}

If D is nonsingular, eliminating y and solving for x yields

\begin{equation} \label{EqSolve.4} {\bf x} = \left( {\bf A} - {\bf B}\,{\bf D}^{-1} {\bf C} \right)^{-1} \left( {\bf p} - {\bf B} \,{\bf D}^{-1} {\bf q} \right) = \left( {\bf M} / {\bf D} \right)^{-1} \left( {\bf p} - {\bf B} \,{\bf D}^{-1} {\bf q} \right) , \end{equation}

where \( {\bf M} / {\bf D} = {\bf A} - {\bf B}\,{\bf D}^{-1} {\bf C} \) is the Schur complement. Similarly, if A is nonsingular, eliminating x and solving for y yields

Solve[{A x + B y == p, CC x + DD y == q}, {x, y}]

{{x -> -((DD p - B q)/(B CC - A DD)), y -> -((-CC p + A q)/(B CC - A DD))}}

\begin{equation} \label{EqSolve.5} {\bf y} = \left( {\bf D} - {\bf C}\,{\bf A}^{-1} {\bf B} \right)^{-1} \left( {\bf q} - {\bf C} \,{\bf A}^{-1} {\bf p} \right) = \left( {\bf M} / {\bf A} \right)^{-1} \left( {\bf q} - {\bf C} \,{\bf A}^{-1} {\bf p} \right) , \end{equation}

where \( {\bf M} / {\bf A} = {\bf D} - {\bf C}\,{\bf A}^{-1} {\bf B} \) is called the Schur complement. An important use of Schur complements is the partitioned solution of linear systems. In fact this is the most important application in Finite Element Method (FEM for short), in connection with superelement analysis.

Example 6: Let us consider a linear algebraic equation M z = r with two input vectors r and b, where

\[ {\bf M} = \begin{bmatrix} 2&5&-2&\phantom{-}1&2 \\ 3&4&-2& \phantom{-}1&2 \\ 1&3&\phantom{-}1&\phantom{-}3&1 \\ 2&1&-3&-2&1 \\ 4&2&-1&\phantom{-}3&1 \end{bmatrix} , \qquad {\bf r} = \begin{bmatrix} 1 \\ 2 \\ 3 \\ 4 \\ 5 \end{bmatrix} \qquad\mbox{and} \qquad {\bf b} = \begin{bmatrix} 1 \\ 7 \\ 3 \\ 4 \\ 5 \end{bmatrix} . \tag{6.1} \]

\[ {\bf M}^{\mathrm T} {\bf y} = {\bf 0} \qquad\mbox{or} \qquad \begin{bmatrix} \phantom{-}2&\phantom{-}3&1&\phantom{-}2&\phantom{-}4 \\ \phantom{-}5&\phantom{-}4&3&\phantom{-}1&\phantom{-}2 \\ -2&-2&1&-3&-1 \\ \phantom{-}1&\phantom{-}1&3&-2&\phantom{-}3 \\ \phantom{-}2&\phantom{-}2&1&\phantom{-}1&\phantom{-}1 \end{bmatrix} \begin{bmatrix} y_1 \\ y_2 \\ y_3 \\ y_4 \\ y_5 \end{bmatrix} = \begin{bmatrix} 0 \\ 0 \\ 0 \\ 0 \\ 0 \end{bmatrix} . \tag{6.2} \]

Equation in blocks

Since the system is too big to handle it manually, we break it into blocks:

\[ {\bf M} = \begin{bmatrix} {\bf A} & {\bf B} \\ {\bf C} & {\bf D} \end{bmatrix} , \tag{6.4} \]

where the partition of the given matrix M into submatrices A, B,C, and D can be chosen in four basic different ways. Of course, there are many other choices because the order of equations does not matter, but we stick with these four choices to demonstrate the method (due to Issai Schur). We consider every case separately.

Four-dimensional matrix A

We consider the case when matrix A is a 4×4 matrix:

\[ {\bf A} = \begin{bmatrix} 2&5&-2&\phantom{-}1 \\ 3&4&-2& \phantom{-}1 \\ 1&3&\phantom{-}1&\phantom{-}3 \\ 2&1&-3&-2 \end{bmatrix} , \qquad {\bf B} = \begin{bmatrix} 2 \\ 2 \\ 1 \\ 1 \end{bmatrix} , \qquad {\bf C} = \begin{bmatrix} 4 & 2 & -1 & 3 \end{bmatrix} , \qquad {\bf D} = \begin{bmatrix} 1 \end{bmatrix} . \tag{6.4.1} \]

Correspondingly, we partition the inbut vectors r and b:

\[ {\bf r} = \begin{bmatrix} {\bf p}_4 \\ {\bf q}_4 \end{bmatrix} , \qquad\mbox{with} \qquad {\bf p}_4 = \begin{bmatrix} 1 \\ 2 \\ 3 \\ 4 \end{bmatrix} , \quad {\bf q}_4 = \begin{bmatrix} 5 \end{bmatrix} , \tag{6.4.2} \]

and

\[ {\bf b} = \begin{bmatrix} {\bf p}_4 \\ {\bf q}_4 \end{bmatrix} , \qquad\mbox{with} \qquad {\bf p}_4 = \begin{bmatrix} 1 \\ 7 \\ 3 \\ 4 \end{bmatrix} , \quad {\bf q}_4 = \begin{bmatrix} 5 \end{bmatrix} . \tag{6.4.3} \]

The system is partitoned as

\[ \begin{cases} 2\,x_1 + 5\,x_2 -2\,x_3 + x_4 + 2\,x_5 &= 1 , \\ 3\,x_1 + 4\,x_2 -2\,x_3 + x_4 + 2\,x_5 &= 2 \quad \mbox{or} \quad 7 , \\ \phantom{-\,}x_1 + 3\,x_2 + x_3 + 3\,x_4 + \phantom{2}x_5 &= 3 , \\ 2\,x_1 + x_2 -3\,x_3 -2\,x_4 +\phantom{2}x_5 &= 4 \\ \hline 4\,x_1 + 2\,x_2 -x_3 + 3\,x_4 + \phantom{2}x_5 &= 5 . \end{cases} \]

This leads to two equations for x and y:

\[ \left( {\bf A} - {\bf B}\,{\bf D}^{-1} {\bf C} \right) {\bf x} = {\bf p} - {\bf B}\,{\bf D}^{-1} {\bf q} \tag{6.4.4} \]

and

\[ \left( {\bf D} - {\bf C}\,{\bf A}^{-1} {\bf B} \right) {\bf y} = {\bf q} - {\bf C}\,{\bf A}^{-1} {\bf p} . \tag{6.4.5} \]

With our partition, these Schur complement matrices become

\[ {\bf M} / {\bf D} = {\bf A} - {\bf B}\,{\bf D}^{-1} {\bf C} = \begin{bmatrix} -12 & -9 & -16 & -13 \\ -11&-10 &-16& -13 \\ -13& -11& -13& -11 \\ -12 & -13&-17&-16 \end{bmatrix} \qquad \Longrightarrow \qquad \left( {\bf M} / {\bf D} \right)^{-1} = \frac{1}{224} \begin{bmatrix} -84& 105& 35& 7 \\ 140&-119& -35& 7 \\ 84& -193& 11& 81 \\ -140&223& 43&-111 \end{bmatrix} \]

A4 = {{2, 5, -2, 1}, {3, 4, -2, 1}, {1, 3, 1, 3}, {2, 1, -3, -2}}
B4 = {2,2,1,1}
C4 = {4,2,-1,3}
MD = A4 - B4.C4

{{-12, -9, -16, -13}, {-11, -10, -16, -13}, {-13, -11, -13, -11}, {-12, -13, -17, -16}}

Its inverse matrix exists because det(M/D) = -224.

Inverse[%]

{{-(3/8), 15/32, -(5/32), 1/32}, {5/8, -(17/32), -(5/32), 1/32}, {3/ 8, -(193/224), 11/224, 81/224}, {-(5/8), 223/224, 43/ 224, -(111/224)}}

but another Schur complement matrix does not exist because matrix A is singular.

Det[A4]

Therefore, our partition does not work: it is impossible to separate the problem into two independent systems of sizes 4 and 1.

Three-dimensional matrix A

We consider the case when matrix A is a 3×3 matrix:

\[ {\bf A} = \begin{bmatrix} 2&5&-2 \\ 3&4&-2 \\ 1&3&\phantom{-}1 \end{bmatrix} , \qquad {\bf B} = \begin{bmatrix} 1 &2 \\ 1&2 \\ 3&1 \end{bmatrix} , \qquad {\bf C} = \begin{bmatrix} 2 & 1 & -3 \\ 4&2&-1 \end{bmatrix} , \qquad {\bf D} = \begin{bmatrix} -2&1 \\ \phantom{-}3&1 \end{bmatrix} . \]

First, we check whether our square matrices are not singular:

A3 = {{2, 5, -2}, {3, 4, -2}, {1, 3, 1}}; Det[A3]

-15

D3 = {{-2, 1}, {3,1}}; Det[D3]

-5

Therefore, it is possible to separate the given 5×5 system of equations into two separate systems of sizes 3×3 and 2×2. The corresponding Schur complements are

\[ {\bf M} / {\bf D} = {\bf A} - {\bf B}\,{\bf D}^{-1} {\bf C} = \begin{bmatrix} -4 & 2 & 2 \\ -3&1 & 2 \\ -3& 1& 2 \end{bmatrix} \]

B3 = {{1, 2}, {1, 2}, {3, 1}}; C3 = {{2,1,-3}, {4,2,-1}};
MD = A3 - B3.Inverse[D3].C3

{{-4, 2, 2}, {-3, 1, 2}, {-3, 1, 2}}

Matrix M/D is singular

Det[MD]

and

\[ {\bf M} / {\bf A} = {\bf D} - {\bf C}\,{\bf A}^{-1} {\bf B} = \frac{1}{3} \begin{bmatrix} 0&\phantom{-}0 \\ 4&-2 \end{bmatrix} . \]

MA = D3 - C3.Inverse[A3].B3

{{0, 0}, {4/3, -(2/3)}}

Matrix M/A is also singular:

Det[MA]

In order to determine solvabilities of each system

\[ \left( {\bf M} / {\bf D} \right) {\bf x} = {\bf p} - {\bf B}\,{\bf D}^{-1} {\bf q} \qquad \mbox{and} \qquad \left( {\bf M} / {\bf A} \right) {\bf y} = {\bf q} -{\bf C}\,{\bf A}^{-1} {\bf p} \tag{6.3.1} \]

we need to find the null-space of each Schur complement.

NullSpace[Transpose[MD]]

and

NullSpace[Transpose[MA]]

So we see that the kernels of Scur complements are spanned on vectors

\[ {\bf y}_3 = \left[ 0, -1, 1 \right]^{\mathrm T} \qquad \mbox{and} \qquad {\bf y}_2 = \left[ 1, 0 \right]^{\mathrm T} , \tag{6.3.2} \]

respectively. Next, we calculate the right-hand sides of Eq.(6.3.1). For vector r, we have

\[ {\bf p} - {\bf B}\,{\bf D}^{-1} {\bf q} = \left[ -8, -7, -2 \right]^{\mathrm T} , \qquad {\bf q} -{\bf C}\,{\bf A}^{-1} {\bf p} = \left[ 5, 1 \right]^{\mathrm T} . \]

p3 = {1, 2, 3} - B3.Inverse[D3].{4, 5}

{-8, -7, -2}

This vector [-8, -7, -2]^T is not orthogonal to y₃ from Eq.(6,3,2). For another 2-dimensional system, we have

{4, 5} - C3.Inverse[A3].{1, 2, 3}

{5, 1}

This vector [5, 1]^T is not orthogonal to y₂ = [1, 0]^T from Eq.(6,3,2), and the system M z = r is inconsistent.

Hoverver, if we choose vector b, situation becomes different.

\[ {\bf p} - {\bf B}\,{\bf D}^{-1} {\bf q} = \left[ -8, -2, -2 \right]^{\mathrm T} , \qquad {\bf q} -{\bf C}\,{\bf A}^{-1} {\bf p} = \left[ 0, -32 \right]^{\mathrm T} . \]

p3 = {1, 7, 3} - B3.Inverse[D3].{4, 5}

{-8, -2, -2}

This vector [-8, -2, -2]^T is orthogonal to y₃ from Eq.(6,3,2). Another two-dimensional vector is

q2 = {4, 5} - C3.Inverse[A3].{1, 7, 3}

{0, -(32/3)}

which is orthogonal to y₂. Therefore, the system M z = b is consistent.

Two-dimensional matrix A

We consider the case when matrix A is chosen as a 2×2 matrix:

\[ {\bf A} = \begin{bmatrix} 2&5 \\ 3&4 \end{bmatrix} , \qquad {\bf B} = \begin{bmatrix} -2&1&2 \\ -2&1&2 \end{bmatrix} , \qquad {\bf C} = \begin{bmatrix} 1 & 3 \\2&1 \\4&2 \end{bmatrix} , \qquad {\bf D} = \begin{bmatrix} \phantom{-}1&\phantom{-}3&1 \\ -3&-2&1 \\ -1& \phantom{-}3&1 \end{bmatrix} . \]

First, we check that our square matrices are not singular:

A2 = {{2, 5}, {3, 4}}; D2 = {{1, 3, 1}, {-3, -2, 1}, {-1, 3, 1}};
Det[D2]

-10

Det[A2]

-7

So we can split the given 5-dimensional systems of equations M z = r and M z = b into two separate equations

\[ \left( {\bf M} / {\bf D} \right) {\bf x} = {\bf p}_2 - {\bf B}\,{\bf D}^{-1} {\bf q}_2 \qquad \mbox{and} \qquad \left( {\bf M} / {\bf A} \right) {\bf y} = {\bf q}_2 -{\bf C}\,{\bf A}^{-1} {\bf p}_2 , \tag{6.2.1} \]

where

\[ {\bf M} / {\bf D} = {\bf A} - {\bf B}\,{\bf D}^{-1} {\bf C} = \begin{bmatrix} -1 & 1 \\ \phantom{-}0&0 \end{bmatrix} \qquad \mbox{and} \qquad {\bf M} / {\bf A} = {\bf D} - {\bf C}\,{\bf A}^{-1} {\bf B} = \frac{1}{7} \begin{bmatrix} \phantom{-}15 & \phantom{-}17& -1 \\ -15& -17 & \phantom{-}1 \\ \phantom{-1}5& \phantom{-}15 & -5 \end{bmatrix} , \tag{6.2.2} \]

B2 = {{-2, 1, 2}, {-2, 1, 2}}; C2 = {{1, 3}, {2, 1}, {4, 2}};
MD = A2 - B2.Inverse[D2].C2

{{-1, 1}, {0, 0}}

MA = D2 - C2.Inverse[A2].B2

{{15/7, 17/7, -(1/7)}, {-(15/7), -(17/7), 1/7}, {5/7, 15/7, -(5/7)}}

and vectors r and b are partitioned as

\[ {\bf r} = \begin{bmatrix} {\bf p} \\ {\bf q} \end{bmatrix} \qquad\mbox{with} \qquad {\bf p}_2 = \begin{bmatrix} 1 \\ 2 \end{bmatrix}, \quad {\bf q}_2 = \begin{bmatrix} 3 \\ 4 \\ 5 \end{bmatrix} \qquad\mbox{and} \qquad {\bf b} = \begin{bmatrix} {\bf p} \\ {\bf q} \end{bmatrix} \qquad\mbox{with} \qquad {\bf p}_2 = \begin{bmatrix} 1 \\ 7 \end{bmatrix}, \quad {\bf q}_2 = \begin{bmatrix} 3 \\ 4 \\ 5 \end{bmatrix} , \]

The right-hand side vectors in Eq.(6.2.1) are

\[ {\bf p}_2 - {\bf B}\,{\bf D}^{-1} {\bf q}_2 = \begin{bmatrix} -6 \\ -5 \end{bmatrix} , \qquad {\bf q}_2 -{\bf C}\,{\bf A}^{-1} {\bf p}_2 = \frac{1}{7} \begin{bmatrix} 18 \\ 17 \\ 13 \end{bmatrix} \tag{6.2.3} \]

{1, 2} - B2.Inverse[D2].{3, 4, 5}

{-6, -5}

{3, 4, 5} - C2.Inverse[A2].{1, 2}

{18/7, 17/7, 13/7}

for vector r and

\[ {\bf p}_2 - {\bf B}\,{\bf D}^{-1} {\bf q}_2 = \begin{bmatrix} -6 \\ \phantom{-}0 \end{bmatrix} , \qquad {\bf q}_2 -{\bf C}\,{\bf A}^{-1} {\bf p}_2 = \frac{1}{7} \begin{bmatrix} \phantom{-}23 \\ -23 \\ -67 \end{bmatrix} \tag{6.2.4} \]

{1, 7} - B2.Inverse[D2].{3, 4, 5}

{-6, 0

{3, 4, 5} - C2.Inverse[A2].{1, 7}

{23/7, -(23/7), -(67/7)}

for input vector b. Finally, we check whether vectors (6.2.3) and (6.2.4) are orthogonal to vectors that span the null space of Schur complements:

\[ {\bf y}_2 = \left[ 0, 1 \right]^{\mathrm T} \qquad \mbox{and} \qquad {\bf y}_3 = \left[ 1, 1, 0 \right]^{\mathrm T} \tag{6.2.5} \]

NullSpace[Transpose[MA]]

NullSpace[Transpose[MD]]

Since vectors (6.2.3), generated by components p and q of the vector r = [p, q]^T, are not orthogonal (their dot products are not zero) to vectors (6.2.5), we claim that system M z = r is inconsistent. On the other hand, vectors (6.2.4), generated by components p = [1, 7]^T and q = ]3, 4, 5]^T of the vector b = [p, q]^T, are orthogonal (their dot products are zeroes) to vectors (6.2.5), and so we claim that system M z = b is consistent.

One-dimensional matrix A

Finally, let us consider the case when matrix A is choisen as a 1×1 matrix:

\[ {\bf A} = \begin{bmatrix} 2 \end{bmatrix} , \qquad {\bf B} = \begin{bmatrix} 5&-2&1&2 \end{bmatrix} , \qquad {\bf C} = \begin{bmatrix} 3 \\ 1 \\ 2 \\ 4 \end{bmatrix} , \qquad {\bf D} = \begin{bmatrix} 4&-2& \phantom{-}1&2 \\ 3&\phantom{-}1&\phantom{-}3&1 \\ 1&-3&-2&1 \\ 2&-1&\phantom{-}3&1 \end{bmatrix} . \]

First, we check whether matrix D is invertible. Mathematica tells us that this is wrong

D1 = {{4, -2, 1, 2}, {3, 1, 3, 1}, {1, -3, -2, 1}, {2, -1, 3, 1}};
Det[D1]

So we conclude that such separation into 1-dimensional and 4-dimensional systems does not lead to an equivalent systems of algebraic equations. ■

R has a dedicated command to solve linear systems of equations \( {\bf A}\,{\bf x} = {\bf b} , \quad \det{\bf A} \ne 0 . \) This means that you can find a solution only for square nonsingular matrix A.

Since the code above requires a package, it will run not directly here (copy onto your own R to run). The following is the output you should expect:

1x1 + 2x2 = 3 

2x1 + 4x2 = 6

x1 + 2x2 = 3 

0 = 0

======================================================

Next, we utilize anothe rpackage plotly for creating interactive web-based graphs via the open source JavaScript:

Rgl is another package that is useful for 3d plotting. It is interactive; however, mac users must download XQuartz to use the rgl package: https://www.xquartz.org

■

Example: We plot a torus with the follwoing commands:

■

 rm(x,y)

https://cran.r-project.org/web/packages/deSolve/vignettes/deSolve.pdf
Soeraert, K., Cash, J., Mazzia, F., Solving differential equations in R, Springer, 2012.
https://www.lehigh.edu/~wes1/R_PDE/R-Matlab_compare.pdf

R TUTORIAL, part 2.1: Linear systems of equations

R solves linear systems