This vignette illustrates the ideas behind solving systems of linear equations of the form Ax=b where
or, spelled out,
[a11a12⋯a1na21a22⋯a2n⋮⋮⋮am1am2⋯amn](x1x2⋮xn)=(b1b2⋮bm) 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)
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.
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.
## 1*x1 - 1*x2 = 2
## 2*x1 + 2*x2 = 1
Check whether they are consistent:
## [1] 2 2
## [1] TRUE
Plot the equations:
## x[1] - 1*x[2] = 2
## 2*x[1] + 2*x[2] = 1
Solve()
is a convenience function that shows the
solution in a more comprehensible form:
## x1 = 5/4
## x2 = -3/4
For three (or more) equations in two unknowns, r(A)≤2, because r(A)≤min. 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.
## 1*x1 - 1*x2 = 2
## 2*x1 + 2*x2 = 1
## 3*x1 + 1*x2 = 3
## [1] 2 2
## [1] TRUE
## x1 = 5/4
## x2 = -3/4
## 0 = 0
Plot the equations:
## x[1] - 1*x[2] = 2
## 2*x[1] + 2*x[2] = 1
## 3*x[1] + x[2] = 3
Three equations in two unknowns are inconsistent when r(\mathbf{A}) < r(\mathbf{A | b}).
## 1*x1 - 1*x2 = 2
## 2*x1 + 2*x2 = 1
## 3*x1 + 1*x2 = 6
## [1] 2 3
## [1] "Mean relative difference: 0.5"
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.
## [,1] [,2] [,3]
## [1,] 1 0 2.75
## [2,] 0 1 -2.25
## [3,] 0 0 -3.00
Solve()
shows this more explicitly, using fractions
where possible:
## x1 = 11/4
## x2 = -9/4
## 0 = -3
An approximate solution is sometimes available using a generalized inverse. This gives \mathbf{x} = (2, -1) as a best close solution.
## [,1]
## [1,] 2
## [2,] -1
Plot the equations. You can see that each pair of equations has a solution, but all three do not have a common, consistent solution.
## x[1] - 1*x[2] = 2
## 2*x[1] + 2*x[2] = 1
## 3*x[1] + x[2] = 6
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.
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)
## 2*x1 + 1*x2 - 1*x3 = 8
## -3*x1 - 1*x2 + 2*x3 = -11
## -2*x1 + 1*x2 + 2*x3 = -3
Are the equations consistent?
## [1] 3 3
## [1] TRUE
Solve for \mathbf{x}.
## x1 x2 x3
## 2 3 -1
Other ways of solving:
## [,1]
## x1 2
## x2 3
## x3 -1
## [,1]
## [1,] 2
## [2,] 3
## [3,] -1
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.
## x1 x2 x3
## [1,] 1 0 0 2
## [2,] 0 1 0 3
## [3,] 0 0 1 -1
`echelon() can be asked to show the steps, as the row operations necessary to reduce \mathbf{X} to the identity matrix \mathbf{I}.
##
## Initial matrix:
## x1 x2 x3
## [1,] 2 1 -1 8
## [2,] -3 -1 2 -11
## [3,] -2 1 2 -3
##
## row: 1
##
## exchange rows 1 and 2
## x1 x2 x3
## [1,] -3 -1 2 -11
## [2,] 2 1 -1 8
## [3,] -2 1 2 -3
##
## multiply row 1 by -1/3
## x1 x2 x3
## [1,] 1 1/3 -2/3 11/3
## [2,] 2 1 -1 8
## [3,] -2 1 2 -3
##
## multiply row 1 by 2 and subtract from row 2
## x1 x2 x3
## [1,] 1 1/3 -2/3 11/3
## [2,] 0 1/3 1/3 2/3
## [3,] -2 1 2 -3
##
## multiply row 1 by 2 and add to row 3
## x1 x2 x3
## [1,] 1 1/3 -2/3 11/3
## [2,] 0 1/3 1/3 2/3
## [3,] 0 5/3 2/3 13/3
##
## row: 2
##
## exchange rows 2 and 3
## x1 x2 x3
## [1,] 1 1/3 -2/3 11/3
## [2,] 0 5/3 2/3 13/3
## [3,] 0 1/3 1/3 2/3
##
## multiply row 2 by 3/5
## x1 x2 x3
## [1,] 1 1/3 -2/3 11/3
## [2,] 0 1 2/5 13/5
## [3,] 0 1/3 1/3 2/3
##
## multiply row 2 by 1/3 and subtract from row 1
## x1 x2 x3
## [1,] 1 0 -4/5 14/5
## [2,] 0 1 2/5 13/5
## [3,] 0 1/3 1/3 2/3
##
## multiply row 2 by 1/3 and subtract from row 3
## x1 x2 x3
## [1,] 1 0 -4/5 14/5
## [2,] 0 1 2/5 13/5
## [3,] 0 0 1/5 -1/5
##
## row: 3
##
## multiply row 3 by 5
## x1 x2 x3
## [1,] 1 0 -4/5 14/5
## [2,] 0 1 2/5 13/5
## [3,] 0 0 1 -1
##
## multiply row 3 by 4/5 and add to row 1
## x1 x2 x3
## [1,] 1 0 0 2
## [2,] 0 1 2/5 13/5
## [3,] 0 0 1 -1
##
## multiply row 3 by 2/5 and subtract from row 2
## x1 x2 x3
## [1,] 1 0 0 2
## [2,] 0 1 0 3
## [3,] 0 0 1 -1
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)
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)
## 1*x1 + 3*x2 + 1*x3 = 2
## 1*x1 - 2*x2 - 2*x3 = 3
## 2*x1 + 1*x2 - 1*x3 = 6
Are the equations consistent? No.
## [1] 2 3
## [1] "Mean relative difference: 0.5"