Example 10: The objective of this example is to demonstrate that No one should have to do an Adjugate or Inverse of a 4 x 4 or hier dimension matrix by hand!!!

We consider a nonsingular 4 × 4 matrix \[ {\bf A} = \begin{bmatrix} -2& \phantom{-}1& -1& \phantom{-}2 \\ -1& \phantom{-}2& -1& \phantom{-}1 \\ \phantom{-}2& -4& \phantom{-}1& -1 \\ \phantom{-}3& \phantom{-}2& \phantom{-}1& -2 \end{bmatrix} . \] Although Mathematica has a dedicated command to calculate the adjugate matrix, we prefer to show all steps in its computation. However, you must be prepared for bewildering, tedious complexity. Calculating the inverse of a 4x4 matrix is a bit laborious but let's break it down into steps. The steps to calculate the inverse of a matrix are as follows:

1. Calculate the determinant of the matrix. If the determinant is 0, the matrix does not have an inverse.
2. Calculate the matrix of minors.
3. Co-factor the matrix of minors.
4. Transpose the co-factor matrix.
5. Divide each element of the transposed matrix by the determinant.

Let's go through each step. Since we have already checked invertibility of matrix A, we go directly to the second step (which is the most time consuming). For a matrix to be invertible it must have a non-zero determinant, so that is the first thing we check

2. Calculate the matrix of minors---it is of size 4 × 4. For each element in the matrix, remove the row and column it's in. The determinant of the resulting 3 × 3 matrix is the minor of that element.

For this exercise to be consistent we need to establish a notation convention. Let's define a symbolic 4 x 4 matrix, M:

Using the convention that the first character is the row number and the second character is the column number, we can double subscript the [lower case] elements of M. Mathematica needs to make these subscripted characters into an object. It does this with the "Symbolize" function in the Notation package. Here is a two-frame animation that cycles the symbolic and numeric repeatedly. Remember the numbers are from our original matrix, A. Notice the pattern of the subscripts as you read left-to-right in a single row; then notice the pattern as you read top-to-bottom in a single column; finally, notice the pattern along the main diagonal from upper left corner to lower right corner.

Having already checked the Determinant, we move to step 2: Calculate the matrix of minors. Each of the 16 cells in the original matrix, A, will have its element replaced by a particular Minor. To find that number, for each element in the matrix, remove the row and column it's in. The determinant of the resulting, smaller, matrix is the minor for that element.

Although Mathematica automatically computes the determinants of each element's minor and places those determinants in a matrix, we perform these operations manually. Once done, you will never want to do it again and you will appreciate use of a computer algebra system.

We need to create a 3x3 "matrix of minors". Imagine if, for each cell, you ignore the other values in the column and row for that cell you would have then remaining a smaller 3x3 matrix: Its determinant is Here is a visual of the minor for the upper left cell in the original matrix. That red square is to be replaced by the Determinant of the minor in the lower right corner. We start with the top left corner and delete first row and first column to generate 3 × 3 matrix

🟥
	2	-1	1
	-4	1	-1
	2	1	-2

The second Minor (Subscript[A, 1,2]) is evaluated as folows. Next, we eliminate the first row and second column from matrix A.

	🟥
-1		-1	1
2		1	-1
3		1	-2

The corresponding matrix \[ {\bf A}_{12} = \begin{bmatrix} -1&-1&\phantom{-}1 \\ \phantom{-}2&\phantom{-}1&-1 \\ \phantom{-}3&\phantom{-}1&-2 \end{bmatrix} , \qquad \det {\bf A}_{12} = -1 , \] has determinant

The next two minors corresponding the first row are \[ {\bf A}_{13} = \begin{bmatrix} -1&\phantom{-}2&-1 \\ \phantom{-}2&-4&-1 \\ \phantom{-}3&\phantom{-}2&-2 \end{bmatrix} , \qquad \det {\bf A}_{13} = 8, \]

		🟥
-1	2		1
2	-4		-1
3	2		-2

and \[ {\bf A}_{14} = \begin{bmatrix} -1&\phantom{-}2&-1 \\ \phantom{-}2&-4&\phantom{-}1 \\ \phantom{-}3&\phantom{-}2&\phantom{-}1 \end{bmatrix} , \qquad \det {\bf A}_{14} = -8 . \]

			🟥
-1	2	-1
2	-4	1
3	2	1

Now we work out the second row.

	1	-1	2
🟥
	-4	1	-1
	2	1	-2

Its determinant is

Minor for (2, 2)-cell is

-2		-1	2
	🟥
2		1	-1
3		1	-2

\[ {\bf A}_{22} = \begin{bmatrix} -2&-1&2 \\ 2&1&-1 \\ 3&1&-2 \end{bmatrix} , \qquad \det {\bf A}_{22} = -1 . \]

Similarly, we have

-2	1		2
		🟥
2	-4		-1
3	2		-2

\[ {\bf A}_{2,3} = \begin{bmatrix} -2 & 1 & 2 \\ 2&-4& -1 \\ 3&2&-2 \end{bmatrix} , \qquad \det {\bf A}_{2,3} = 13 , \]

-2	1	-1
			🟥
2	-4	1
3	2	1

\[ {\bf A}_{2,4} = \begin{bmatrix} -2 & 1 & -1 \\ 2&-4& 1 \\ 3&2&1 \end{bmatrix} , \qquad \det {\bf A}_{2,4} = -3 , \]

Reaching the halfway point, we can again check against Mathematica's built in function at this stage. First, let's assemble our individual minors.

There is some consistency, albeit little. Notice the 2 x 4 matrix above is the reverse and reordered of the last two lines in the 4 x 4 matrix below We can force them to match

The ninth Minor - third row, first column is

	1	-1	2
	2	-1	1
🟥
	2	1	-2

\[ {\bf A}_{31} = \begin{bmatrix} 1&-1&\phantom{-}2 \\ 2&-1&\phantom{-}1 \\ 2& \phantom{-}1&-2 \end{bmatrix} , \qquad \det {\bf A}_{31} = 3, \]

-2		-1	2
-1		-1	1
	🟥
3		1	-2

\[ {\bf A}_{32} = \begin{bmatrix} -2& -1& \phantom{-}2 \\ -1& -1& \phantom{-}1 \\ \phantom{-}3& \phantom{-}1& -2 \end{bmatrix} , \qquad \det {\bf A}_{32} = 1, \]

-2	1		2
-1	2		1
		🟥
3	2		-2

\[ {\bf A}_{33} = \begin{bmatrix} -2&1&\phantom{-}2 \\ -1&2&\phantom{-}1 \\ \phantom{-}3&2&-2 \end{bmatrix} , \qquad \det {\bf A}_{33} = -3 , \]

-2	1	-1
-1	2	-1
			🟥
3	2	1

The twelfth Minor - third row, fourth column is \[ {\bf A}_{34} = \begin{bmatrix} -2&1&-1 \\ -1&2&\phantom{-}1 \\ \phantom{-}3&2&\phantom{-}1 \end{bmatrix} , \qquad \det {\bf A}_{34} = -2 , \] The thirteenth Minor - fourth row, first column (Subscript[A, 4,1]) is

	1	-1	2
	-4	1	-1
	2	1	-2
🟥

\[ {\bf A}_{41} = \begin{bmatrix} \phantom{-}1&-1&\phantom{-}2 \\ \phantom{-}2&-1&\phantom{-}1 \\ -4&\phantom{-}1&-1 \end{bmatrix} , \qquad \det {\bf A}_{41} = -2 , \]

-2		-1	2
-1		-1	1
2		1	-1
	🟥

\[ {\bf A}_{42} = \begin{bmatrix} -2&-1&\phantom{-}2 \\ -1&-1&\phantom{-}1 \\ \phantom{-}2&\phantom{-}1&-1 \end{bmatrix} , \qquad \det {\bf A}_{42} = 1 , \] The fifteenth Minor - fourth row, third column (Subscript[A, 4,3]) is

-2	1		2
-1	2		1
2	-4		-1
		🟥

\[ {\bf A}_{43} = \begin{bmatrix} -2&\phantom{-}1&\phantom{-}2 \\ -1&\phantom{-}2&\phantom{-}1 \\ \phantom{-}2&-4&-1 \end{bmatrix} , \qquad \det {\bf A}_{43} = -3 , \] The sixteenth Minor - fourth row, fourth column (Subscript[A, 4,4]) is

-2	1	-1
-1	2	-1
2	-4	1
			🟥

\[ {\bf A}_{44} = \begin{bmatrix} -2&\phantom{-}1&-1 \\ -1&\phantom{-}2&-1 \\ \phantom{-}2&-4&\phantom{-}1 \end{bmatrix} , \qquad \det {\bf A}_{44} = 3 . \]

Finally, we are now in a position to make the matrix of minors for our original matrix A based on the names we have given each one. We than can compare it with the built-in Mathematica code for Minors[ ]:

\[ {\bf M} = \begin{bmatrix} \phantom{-}2&-1& \phantom{-}8&-8 \\ -3&-1&13 &-3 \\ \phantom{-}3&\phantom{-}1&-3&-2 \\ -2&\phantom{-}1&-3&\phantom{-}3 \end{bmatrix} . \]

With some machinations we can force Mathematica to agree with us!

Or, we can achieve the same answer by contorting our matrix of minors (minMat) equivalently.

3. Co-factor the matrix of minors. Apply a checkerboard of minuses to the "Matrix of Minors". Starting with a plus for the top left, alternate plus and minus through the matrix.

\[ \mbox{Co-factors} = \begin{bmatrix} \phantom{-}2&1& \phantom{-}8&8 \\ 3&-1&-13 &-3 \\ \phantom{-}3&-1&-3&2 \\ 2&\phantom{-}1&3&\phantom{-}3 \end{bmatrix} . \]

4. Transpose the co-factor matrix. Swap the elements of the rows and columns.

\[ \mbox{Co-factors}^{\mathrm T} = \begin{bmatrix} 2&\phantom{-}3& \phantom{-}3&2 \\ 1&-1&-1 &1 \\ 8&-13&-3&3 \\ 8& -3&\phantom{-}2&3 \end{bmatrix} . \]

5. Divide each element of the transposed matrix by the determinant. This will result in the inverse of the original matrix.

We compare our answer with the answer provided by Mathematica: \[ \mbox{Adjugate} = \left( \det {\bf A} \right) {\bf A}^{-1} = \begin{bmatrix} 2& \phantom{-}3& \phantom{-}3& 2 \\ 1& -1& -1& 1 \\ 8& -13& -3& 3 \\ 8& -3& \phantom{-}2& 3\end{bmatrix} . \]

Example 11: We start with a general 2 × 2 matrix: and where we use the vertical lines to define determinants of corresponding matrices.

The cofactor matrix can be computed in the Wolfram language using

which is the equivalent of the (i,j)-th components of the CofactorMatrix defined below:

Example 12: Consider the matrix To find its inverse, we need to check whether the given matrix is singular or not: Now we calculate the cofactor matrix: Similarly, we can calculate all other cofactors, which yields, together with \( \det{\bf A} = -1 , \) the inverse matrix:

Example 13: Suppose we need to find the inverse of the matrix Of course, Mathematica finds its inverse in a blank of eye: \[ {\bf A}^{-1} = \frac{1}{2} \begin{bmatrix} \phantom{-}3&-1&6 \\ -7&\phantom{-}3&-20 \\ \phantom{-}3&-1&8 \end{bmatrix} . \tag{13.1} \] We determine every column in the latter by solving three systems of equations: \[ \begin{split} 2\, x_1 + x_2 + x_3 &= 1, \\ -2\, x_1 + 3\, x_2 + 9\, x_3 &= 0 , \\ -x_1 + x_3 &= 0 ; \end{split} \tag{13.2} \] According to the theory, to solve system of equations (13.2), we build the augmented matrix \[ \left[ {\bf A} \mid {\bf b} \right] = \begin{bmatrix} \phantom{-}2 & 1 & 1&1 \\ -2& 3 &9&0 \\ -1&0& 1&0 \end{bmatrix} . \] Then we use elementary row operations to reduce the augmented matrix to upper traingular form. \[ \left[ {\bf A} \mid {\bf b} \right] \,\underset{R}{\sim}\,\begin{bmatrix} 2 & 1 & 1&1 \\ 0& 4 &10&1 \\ 0&\frac{1}{2}& \frac{3}{2}&\frac{1}{2} \end{bmatrix} \,\underset{R}{\sim}\,\begin{bmatrix} 2 & 1 & 1&1 \\ 0& 4 &10&1 \\ 0&0& \frac{1}{4}&\frac{3}{8} \end{bmatrix} . \] Therefore, system of equations (13.2) is equivalent to \[ \begin{split} 2\, x_1 + x_2 + x_3 &= 1, \\ 4\, x_2 + 10\, x_3 &= 1 , \\ \frac{1}{4}\, x_3 &= \frac{3}{8} . \end{split} \] From the last equation, we immediately read off x₃ = 3/2. Substituting this value into the rest two equations, we get \[ \begin{split} 2\, x_1 + x_2 + \frac{3}{2} &= 1, \\ 4\, x_2 + 10\, \frac{3}{2} &= 1 , \end{split} \qquad \Longrightarrow \qquad \begin{split} 2\, x_1 + x_2 &= -\frac{1}{2}, \\ 4\, x_2 &= -14 . \end{split} \] Hence, the second unknown becomes x₂ = −7/2, which we use to determine the last one \[ 2\, x_1 - \frac{7}{2} = - \frac{1}{2} \qquad \Longrightarrow \qquad x_1 = \frac{3}{2} . \] Finally, we obtain the first column of the inverse matrix A⁻¹ to be \[ {\bf c}_1 = \frac{1}{2}\,\begin{bmatrix} 3 \\ -7 \\ 3 \end{bmatrix} . \] Then we come to the second system of equations: \[ \begin{split} 2\, x_1 + x_2 + x_3 &= 0, \\ -2\, x_1 + 3\, x_2 + 9\, x_3 &= 1 , \\ -x_1 + x_3 &= 0 ; \end{split} \tag{13.3} \] Determination of its solution we deligate to Mathematica: So the second column of the inverse matrix becomes \[ {\bf c}_2 = \frac{1}{2}\,\begin{bmatrix} -1 \\ 3 \\ -1 \end{bmatrix} . \] We similarly solve the system \[ \begin{split} 2\, x_1 + x_2 + x_3 &= 0, \\ -2\, x_1 + 3\, x_2 + 9\, x_3 &= 0 , \\ -x_1 + x_3 &= 1 ; \end{split} \tag{13.4} \] Hence, the last column becomes \[ {\bf c}_3 =\begin{bmatrix} 3 \\ -10 \\ 4 \end{bmatrix} . \] Finally, we build the inverse matrix from these column vectors: \[ {\bf A}^{-1} = \left[ {\bf c}_1 \ {\bf c}_2 \ {\bf c}_3 \right] . \]

Example 14: Let us consider the defective 4 × 4 matrix \[ {\bf B} = \begin{bmatrix} -115& 17& -85& -15 \\ -62& 11& -45& -8 \\ 124& -18& 92& 16 \\ 130& -19& 95& 18 \end{bmatrix} . \tag{14.1} \] First, we find its resolvent using Mathematica \[ {\bf R}_{\lambda}({\bf B}) = \left( \lambda {\bf I} - {\bf B} \right)^{-1} = \frac{m(\lambda )}{\left( \lambda -1 \right)^2 \left( \lambda -2 \right)} , \] where the numerator matrix in the resolvent formula is \[ m(\lambda ) = \begin{bmatrix} 146 - 119 \lambda + \lambda^2 & -21 + 17 \lambda & 105 - 85 \lambda & 19 - 15 \lambda \\ 76 - 62 \lambda & -10 + 7 \lambda + \lambda^2 & 55 - 45 \lambda & 10 - 8 \lambda \\ -152 + 124 \lambda & 22 - 18 \lambda &-109 + 88 \lambda + \lambda^2 & -20 + 16 \lambda \\ -172 + 130 \lambda & 25 - 19 \lambda & -125 + 95 \lambda & -21 + 14 \lambda + \lambda^2 \end{bmatrix} . \] Since the resolvent has two singular points λ = 1 and λ = 2, we need to find two residues at these points. We start with the first one. \[ \mbox{Res}_{\lambda =1} \frac{m(\lambda )}{\lambda \left( \lambda -1 \right)^2 \left( \lambda -2 \right)} = \lim_{\lambda \to 1} \frac{\text d}{{\text d}\lambda} \left( \frac{m(\lambda )}{\lambda \left( \lambda -2 \right)} \right) . \] This limit can be evaluated with the aid of Mathematica: So the residue at λ = 1 becomes \[ {\bf P} = \begin{bmatrix} 117& -17& 85& 15 \\ 62& -9& 45& 8 \\ -124& 18& -90& -16 \\ -130& 19& -95& -16 \end{bmatrix} . \] The residue at another singular point is \[ \mbox{Res}_{\lambda =2} \frac{m(\lambda )}{\lambda \left( \lambda -1 \right)^2 \left( \lambda -2 \right)} = \lim_{\lambda \to 2} \frac{m(\lambda )}{\lambda \left( \lambda -1 \right)^2} = \frac{m(2)}{2} . \] Of course, this value is not hard to evaluate by hand, but lazy people like me deligate this job to software: \[ \mbox{Res}_{\lambda =2} \frac{m(\lambda )}{\lambda \left( \lambda -1 \right)^2 \left( \lambda -2 \right)} = \begin{bmatrix} -44& 13/2& -(65/2)& -(11/2) \\ -24& 4& -(35/2)& -3 \\ 48& -7& 71/2& 6 \\ 44& -(13/2)& 65/2& 11/2 \end{bmatrix} . \] Adding these two expressions, we obtain the explicit expression for the inverse matrix: We compare this answer with the output provided by Mathematica when it was asked to find B⁻¹: \[ {\bf B}^{-1} = \begin{bmatrix} 73 & - \frac{21}{2} & \frac{105}{2} & \frac{19}{2} \\ 38 & -5 & \frac{55}{2} & 5 \\ -76 & 11 & -\frac{109}{2} &-10 \\ -86 & \frac{25}{2} & - \frac{125}{2} &- \frac{21}{2} \end{bmatrix} . \]

Example 15: Consider the spring-mass system that contains two masses m₁ and m₂ connecting with each other by a spring and also connected withe wall by two prings. This system is depicted below. Let x₁ and x₂ be displacements from the equilibrium positions of masses m₁ and m₂, respectively (with positive displacement being to the right). We assume that these masses are powered by forces f₁ and f₂, pointed positive when they act to the right.

According to Hooke’s Law, the force F required to stretch or compress a spring with spring constant k by a distance x is given by \[ F = k\,x . \] Taking into account the direction of force applied to each mass by the spring, we can write an equation for the force acting on each spring. The mass m₁ is being acted upon by two springs with spring constants k₁ and k₂. A displacement of x₁ by mass m₁ will cause the spring k₁ to stretch or compress depending on if x₁ is positive or negative. In either case, by Hooke’s Law, the force acting upon mass m₁ by spring k₁ is \[ - k_1 x_2 . \] We needed to include the negative sign as if x₁ were positive, then the force will need to be negative as the spring will be trying to pull the mass back.

How much spring k₂ will be stretched or compressed depends on displacements x₁ and x₂. In particular, it will be stretched by x₂ − x₂. Hence, the force acting on mass m₁ by spring k₂ is \[ k_2 \left( x_2 - x_1 \right) . \] For the system to be in equilibrium, the sum of the forces acting on mass m₁ must equal 0. Thus, we have \[ f_1 - k_1 x_1 + k_2 \left( x_2 - x_1 \right) = 0 . \]

Mass m₂ is being acted upon by springs k₂ and k₃. Using similar analysis, we find that \[ f_2 - k_3 x_2 -k_2 \left( x_2 - x_1 \right) = 0 . \] Rearranging this equations gives \[ \begin{split} k_1 x_1 + k_2 \left( x_1 - x_2 \right) &= f_1 , \\ k_2 \left( x_2 - x_1 \right) + k_3 x_2 &= f_1 . \end{split} \]

We want to find the forces f₁ and f₂ acting on the masses m₁ and m₂ given displacements x₁ and x₂. By Hooke’s Law, we find that the forces f₁ and f₂ satisfy Using matrix/vector multiplication, we can rewrite this system as \[ {\bf K}\,{\bf x} = {\bf f}, \qquad {\bf K} = \begin{bmatrix} k_1 + k_2 & - k_2 \\ -k_2 & k_2 + k_3 \end{bmatrix}, \quad {\bf x} = \begin{bmatrix} x_1 \\ x_2 \end{bmatrix}, \quad {\bf f} = \begin{bmatrix} f_1 \\ f_2 \end{bmatrix} . \] Matrix K is known as the stiffness matrix of the system. In order to find displacements, we multiply both sides of the latter equation by K⁻¹ to obtain \[ {\bf K}^{-1} = \frac{1}{k_1 k_2 + k_1 k_3 + k_2 k_3} \begin{bmatrix} k_2 + k_3 & k_2 \\ k_2 & k_1 + k_2 \end{bmatrix} . \]

1. Anton, Howard (2005), Elementary Linear Algebra (Applications Version) (9th ed.), Wiley International
2. Beezer, R.A., A First Course in Linear Algebra, 2017.
3. Dobrushkin, V.A., Applied Differential Equations. The Primary Course, second edition, CRC Press2022.
4. Fadeev--LeVerrier algorithm, Wikipedia.
5. Frame, J.S., A simple recursion formula for inverting a matrix, Bulletin of the American Mathematical Society, 1949, Vol. 55, p. 1045. doi:10.1090/S0002-9904-1949-09310-2
6. Greenspan, D., Methods of matrix inversion, The American mathematical Monthly, 1955, Vol. 62, No. pp. 303--318.
7. Karlsson, L., Computing explicit matrix inverses by recursion, Master's Thesis in Computing Science, 2006.
8. Lightstone, A.H., Two methods of inverting matrices, Mathematics Magazine, 1968, Vol. 41, No. 1, pp. 1--7.