knitr::opts_chunk$set(
  warning = FALSE,
  message = FALSE
)
options(digits=4)
par(mar=c(5,4,1,1)+.1)
library(rgl)
library(knitr)
knit_hooks$set(webgl = hook_webgl)
library(matlib)   # use the package

This vignette illustrates the ideas behind solving systems of linear equations of the form $\mathbf{A x = b}$ where

or, spelled out,

$$ \begin{bmatrix} a_{11} & a_{12} & \cdots & a_{1n} \ a_{21} & a_{22} & \cdots & a_{2n} \ \vdots & \vdots & & \vdots \ a_{m1} & a_{m2} & \cdots & a_{mn} \ \end{bmatrix} \begin{pmatrix} x_{1} \ x_{2} \ \vdots \ x_{n} \ \end{pmatrix} \quad=\quad \begin{pmatrix} b_{1} \ b_{2} \ \vdots \ b_{m} \ \end{pmatrix} $$ For three equations in three unknowns, the equations look like this:

A <- matrix(paste0("a_{", outer(1:3, 1:3, FUN  = paste0), "}"), 
            nrow=3) 
b <- paste0("b_", 1:3)
x <- paste0("x", 1:3)
showEqn(A, b, vars = x, latex=TRUE)

Conditions for a solution

The general conditions for solutions are:

We use c( R(A), R(cbind(A,b)) ) to show the ranks, and all.equal( R(A), R(cbind(A,b)) ) to test for consistency.

library(matlib)   # use the package

Equations in two unknowns

Each equation in two unknowns corresponds to a line in 2D space. The equations have a unique solution if all lines intersect in a point.

Two consistent equations

A <- matrix(c(1, 2, -1, 2), 2, 2)
b <- c(2,1)
showEqn(A, b)

Check whether they are consistent:

c( R(A), R(cbind(A,b)) )          # show ranks
all.equal( R(A), R(cbind(A,b)) )  # consistent?

Plot the equations:

#| fig.alt: Plot of two consistent equations which plot as lines intersecting in a point
par(mar=c(4,4,0,0)+.1)
plotEqn(A,b)

Solve() is a convenience function that shows the solution in a more comprehensible form:

Solve(A, b, fractions = TRUE)

Three consistent equations

For three (or more) equations in two unknowns, $r(\mathbf{A}) \le 2$, because $r(\mathbf{A}) \le \min(m,n)$. The equations will be consistent if $r(\mathbf{A}) = r(\mathbf{A | b})$. This means that whatever linear relations exist among the rows of $\mathbf{A}$ are the same as those among the elements of $\mathbf{b}$.

Geometrically, this means that all three lines intersect in a point.

A <- matrix(c(1,2,3, -1, 2, 1), 3, 2)
b <- c(2,1,3)
showEqn(A, b)
c( R(A), R(cbind(A,b)) )          # show ranks
all.equal( R(A), R(cbind(A,b)) )  # consistent?

Solve(A, b, fractions=TRUE)       # show solution 

Plot the equations:

#| fig.alt: Plot of three consistent equations which plot as three lines intersecting in a point
par(mar=c(4,4,0,0)+.1)
plotEqn(A,b)

Three inconsistent equations

Three equations in two unknowns are inconsistent when $r(\mathbf{A}) < r(\mathbf{A | b})$.

A <- matrix(c(1,2,3, -1, 2, 1), 3, 2)
b <- c(2,1,6)
showEqn(A, b)
c( R(A), R(cbind(A,b)) )          # show ranks
all.equal( R(A), R(cbind(A,b)) )  # consistent?

You can see this in the result of reducing $\mathbf{A} | \mathbf{b}$ to echelon form, where the last row indicates the inconsistency because it represents the equation $0 x_1 + 0 x_2 = -3$.

echelon(A, b)

Solve() shows this more explicitly, using fractions where possible:

Solve(A, b, fractions=TRUE)

An approximate solution is sometimes available using a generalized inverse. This gives $\mathbf{x} = (2, -1)$ as a best close solution.

x <- MASS::ginv(A) %*% b
x

Plot the equations. You can see that each pair of equations has a solution, but all three do not have a common, consistent solution.

#| fig.alt: Plot of the lines corresponding to three inconsistent equations. They do not all intersect in a point, indicating that there is no common solution.
par(mar=c(4,4,0,0)+.1)
plotEqn(A,b, xlim=c(-2, 4))
# add the ginv() solution
points(x[1], x[2], pch=15)

Equations in three unknowns

Each equation in three unknowns corresponds to a plane in 3D space. The equations have a unique solution if all planes intersect in a point.

Three consistent equations

An example:

A <- matrix(c(2, 1, -1,
             -3, -1, 2,
             -2,  1, 2), 3, 3, byrow=TRUE)
colnames(A) <- paste0('x', 1:3)
b <- c(8, -11, -3)
showEqn(A, b)

Are the equations consistent?

c( R(A), R(cbind(A,b)) )          # show ranks
all.equal( R(A), R(cbind(A,b)) )  # consistent?

Solve for $\mathbf{x}$.

solve(A, b)

Other ways of solving:

solve(A) %*% b
inv(A) %*% b

Yet another way to see the solution is to reduce $\mathbf{A | b}$ to echelon form. The result of this is the matrix $[\mathbf{I \quad | \quad A^{-1}b}]$, with the solution in the last column.

echelon(A, b)

`echelon() can be asked to show the steps, as the row operations necessary to reduce $\mathbf{X}$ to the identity matrix $\mathbf{I}$.

echelon(A, b, verbose=TRUE, fractions=TRUE)

Now, let's plot them.

plotEqn3d() uses rgl for 3D graphics. If you rotate the figure, you'll see an orientation where all three planes intersect at the solution point, $\mathbf{x} = (2, 3, -1)$

plotEqn3d(A,b, xlim=c(0,4), ylim=c(0,4))

Three inconsistent equations

A <- matrix(c(1,  3, 1,
              1, -2, -2,
              2,  1, -1), 3, 3, byrow=TRUE)
colnames(A) <- paste0('x', 1:3)
b <- c(2, 3, 6)
showEqn(A, b)

Are the equations consistent? No.

c( R(A), R(cbind(A,b)) )          # show ranks
all.equal( R(A), R(cbind(A,b)) )  # consistent?


friendly/matlib documentation built on Dec. 6, 2024, 9:25 a.m.