R has sufficient built-in facilities to solve systems of linear equations numerically. We illustrate them using the system \[ \left\{\begin{array}{rrrrrl} x & + & 2y & + & 3z & =\;2\\ & & y & + & 2z & =\;1\\ 3x & + & y & + & z & =\;3 \end{array}\right. \] with one solution, namely \[ \cv{x\\y\\z}=
\left(\!\!\begin{array}{r} 1\\-1\\1\end{array}\right)\]

We start with entering the system of equations in an R session in the appropriate matrix form.

> A <- matrix(c(1, 2, 3,
+               0, 1, 2,
+               3, 1, 1), nrow=3, byrow=TRUE)
> A
     [,1] [,2] [,3]
[1,]    1    2    3
[2,]    0    1    2
[3,]    3    1    1
> b <- matrix(c(2, 1, 3), nrow=3); b
     [,1]
[1,]    2
[2,]    1
[3,]    3

We can now use solve to find the vector \(\vec{v}\) so that \(A\vec{v}=\vec{b}\)

> v <- solve(A,b); v
     [,1]
[1,]    1
[2,]   -1
[3,]    1

Advantage of the matrix form is that you can better control the theory for solving equations. For example, you can apply Gaussian elimination to solve the system. The name of the function rref from the pracma package is an acronym for reduced row echelon form.

> Ab <- cbind(A,b);  Ab  # de gerande matrix
     [,1] [,2] [,3] [,4]
[1,]    1    2    3    2
[2,]    0    1    2    1
[3,]    3    1    1    3
> library(pracma)
> M <- rref(Ab); M  # berekening van de gereduceerde trapvorm
     [,1] [,2] [,3] [,4]
[1,]    1    0    0    1
[2,]    0    1    0   -1
[3,]    0    0    1    1
> M(,4)  # selectie van de vierde kolom
[1] 1 -1 1