In this article, we discus systems of linear equations as treated in linear algebra.

```
import numpy as np
import matplotlib.pyplot as plt
from sympy import *
from sympy.plotting import *
x, y, z = symbols("x y z")
```

## Linear Equations in \(n\) Variables

In analytical geometry, we learn that a line in two-dimensional space has the standard form

\[a_1x + a_2y = b\]

\(a_1\), \(a_2\), and \(b\) are constants. This line is a linear equation in two variables (hence why it is in two-dimensional space). This fact extends to three-dimensional spaces as well.

\[a_1x + a_2y + a_3z = b\]

By definition, then, a **linear equation in \(n\) variables** \(x_1, x_2, x_3, \cdots, x_n\) has the form

\[a_1x_1 + a_2x_2 + a_3x_3 + \cdots + a_nx_n = b\]

The **coefficients** \(a_1, a_2, \cdots, a_n\) are real numbers, and the **constant term** \(b\) is a real number. \(a_1\) and \(x_1\) are, respectively, the **leading coefficient** and the **leading variable**. Linear equations have variables that only appear to the first power. Their variables are not involved in trigonometric, exponential, or logarithmic functions.

A **solution** of a linear equation in \(n\) variables is a sequence of \(n\) real numbers \(s_1, s_2, \cdots, s_n\) so that the equation is true when they’re substituted for the variables of the equation.

Let’s solve the linear equation \(3x + 2y - z = 3\) for each variable.

```
eq = Eq(3*x + 2*y - z, 4)
solutions = [solve(eq, x, dict=True), solve(eq, y, dict=True), solve(eq, z, dict=True)]
solutions
```

`## [[{x: -2*y/3 + z/3 + 4/3}], [{y: -3*x/2 + z/2 + 2}], [{z: 3*x + 2*y - 4}]]`

### Systems of Linear Equations

A **system of \(m\) linear equations in \(n\) variables** is a set of \(m\) equations, each having the same \(n\) variables.

\[ a_{11}x_1 + a_{12}x_2 + a_{13}x_3 + \cdots + a_{1n}x_n = b_1 \\ a_{21}x_1 + a_{22}x_2 + a_{23}x_3 + \cdots + a_{2n}x_n = b_2 \\ \vdots \\ a_{m1}x_1 + a_{m2}x_2 + a_{m3}x_3 + \cdots + a_{mn}x_n = b_n \\ \]

A **solution** of a system of linear equations is a sequence of numbers \(s_1, s_2, \cdots, s_n\) that is a solution of each of the linear equations in the system. This solution, geometrically and conventionally, corresponds to the point where the lines of equations meet.

A system of linear equations can have one solution (if the lines intersect at a single point), in which case it’s called a **consistent system**, infinite solutions (if the lines are “on top of each other”, i.e. they are coincident lines), or no solution (if they are parallel), in which case it’s called an **inconsistent system**.

### Systems of Linear Equations in Two Variables

Let’s graph and solve the following linear equation.

\[ x + y = 3 \\ x - y = -1 \]

```
# rewrite functions by moving the right hand side to the left hand side
equations = [x + y - 3, x - y + 1]
linsolve(equations, x, y)
```

`## FiniteSet((1, 2))`

```
# solve each function to get a form we can plot with sympy
a = solve(equations[0], x)
a
```

`## [3 - y]`

```
b = solve(equations[1], x)
b
```

`## [y - 1]`

```
p = plot(a[0], line_color='b', legend=True, show=False)
p.extend(plot(b[0], line_color='r', legend=True, show=False))
p.show()
```

## Solving Systems of Linear Equations

Equations in Row-Echelon form (meaning, written in a “triangular”, stair-step pattern with leading coefficients of 1) can be solved using Back-Substitution. Systems of equations that are not in Row-Echelon form can be re-written into that form using Gaussian Eliminatoin. Two systems of linear equations are **equivalent** if they have the same solution set.

Let’s solve the system

\[ x - 2y + 3z = 9 \\ -x + 3y = -4 \\ 2x - 5y + 5z = 17 \]

```
equations = [x - 2*y + 3*z -9, -x + 3*y + 4, 2*x -5*y + 5*z -17]
linsolve(equations, x, y, z)
```

`## FiniteSet((1, -1, 2))`

## Introduction to Matrices

A **matrix** is an array of objects (numbers, expressions, or symbols) arranged in a rectangular shape. We denote matrices are \(m \times n\) matrices, where \(m\) denotes the number of rows and \(n\) denotes the number of columns. \(m\) and \(n\) are what determines the **size** of the matrix. We denote each entry of the matrix as \(a_{ij}\) which indicates its row index (\(i\)) and column index (\(j\)).

One common use of matrices is to represent systems of linear equations. The matrix derived from the coefficients and constant terms of a system is called the **augmented matrix** of the system. The matrix containing only the coefficients is the **coefficient matrix**. The methods of Back-Substitution and Gaussian / Gauss-Jordan Elimination can be used on matrices as well. The elimination operations turn a matrix into Row-Echelon form for later Back-Substitution. Although, we can use programs to solve systems of linear equations, let’s illustrate the Row-Echelon form.

Let’s put the system of equations form above into Row-Echelon form.

\[ x - 2y + 3z = 9 \\ -x + 3y = -4 \\ 2x - 5y + 5z = 17 \]

The augmented matrix representation of this system is

\[ \begin{bmatrix} \phantom{-}1 & \phantom{-}2 & 3 & \phantom{-}9 \\ -1 & \phantom{-}3 & 0 & -4 \\ \phantom{-}2 & -5 & 5 & \phantom{-}17 \end{bmatrix} \]

```
## variables with no coefficient get a zero.
A = Matrix([[1, -2, 3, 9], [1, 3, 0, -4], [2, -5, 5, 17]])
A
```

```
## Matrix([
## [1, -2, 3, 9],
## [1, 3, 0, -4],
## [2, -5, 5, 17]])
```

`A.rref()[0]`

```
## Matrix([
## [1, 0, 0, -1/4],
## [0, 1, 0, -5/4],
## [0, 0, 1, 9/4]])
```

## Applications of Systems of Linear Equations

Systems of linear equations have a variety of applications, we’re treating a few ones below.

### Polynomial Curve Fitting

Suppose we have a collection of data represented by \(n\) points in the \(xy\)-plane,

\[(x_1, y_1), (x_2, y_2) \cdots, (x_n, y_n)\]

and we want to find a polynomial function with a degree \(n - 1\)

\[p(x) = a_0 + a_1x + a_2x^2 + \cdots + a_{n - 1}x^{n - 1}\]

whose graph passes through the points. This is what we called **polynomial curve fitting**. To solve for the \(n\) coefficients of \(p(x)\), we substitute each of the \(n\) points into the polynomial function and obtain \(n\) linear equations in \(n\) variables \(a_0, a_1, a_2, \cdots, a_{n-1}\).

As an example, let’s determine the polynomial \(p(x) = a_0 + a_1x + a_2x^2\) that passes through the points \((1, 4), (2, 0), (3, 12)\). Using the aforementioned substitution we get

\[ p(1) = a_0 + a_1(1) + a_2(1)^2 = a_0 + a_1 + a_2 = 4 \\ p(2) = a_0 + a_1(2) + a_2(2)^2 = a_0 + 2a_1 + 4a_2 = 0 \\ p(3) = a_0 + a_1(3) + a_2(3)^2 = a_0 + 3a_1 + 9a_2 = 12 \]

```
A = Matrix([[1, 1, 1, 4], [1, 2, 4, 0], [1, 3, 9, 12]])
solve_linear_system_LU(A, [x, y, z])
```

`## {x: 24, y: -28, z: 8}`

The solution that tells us that the polynomial function is

\[p(x) = 24 - 28x + 8x^2\]

```
def p(x):
return 24 - 28*x + 8*x**2
x = np.linspace(0, 4)
y = p(x)
plt.plot([1, 2, 3], [4, 0, 12], 'ro')
plt.plot(x, y)
```

Let’s illustrate an application from astrophysics. Let’s find the polynomial that relates the periods of the first three planets to their mean distances from the sun. Then, let’s test the accuracy of the fit by using the polynomial to calculate the period of Mars.

Mercury has a mean distance of 0.387 and period of 0.241. Venus has a mean distance of 0.723 and period of 0.615. Earth has a mean distance of 1.0 and period of 1.0.

Our polynomial will have to be of degree \(n - 1 = 3 - 1 = 2\).

\[p(x) = a_0 + a_1x + a_2x^2\]

fitted to the points \((0.387, 0.241), (0.723, 0.615), (1, 1)\). This gives us a system of liner equation

\[ a_0 + 0.387a_1 + (0.387)^2a_2 = 0.241 \\ a_0 + 0.723a_1 + (0.723)^2a_2 = 0.615 \\ a_0 + a_1 + a_2 = 1 \]

```
a0, a1, a2 = symbols("a0 a1 a2")
solve_linear_system_LU(Matrix([
[1, 0.387, 0.387**2, 0.241],
[1, 0.723, 0.723**2, 0.615],
[1, 1, 1, 1]]), [a0, a1, a2])
```

`## {a0: -0.0634254004898172, a1: 0.611881422258717, a2: 0.451543978231100}`

Therefore the function is

\[p(x) = -0.0634 + 0.6119x + 0.4515x^2\]

To find a polynomial fit without constructing the whole matrix, we can use `numpy`

and encode the \(xy\)-coordinates into arrays. The `np.polyfit`

function returns the coefficients with the highest power first.

```
x = np.array([0.387, 0.723, 1])
y = np.array([0.241, 0.615, 1])
np.polyfit(x, y, 2)
```

`## array([ 0.45154398, 0.61188142, -0.0634254 ])`

Knowing that mars has a mean distance of 1.523 from the sun

```
def p(x):
return -0.0634 + 0.6119*x + 0.4515*x**2
# mars period
p(1.523)
```

`## 1.9157910434999998`

```
from sympy import *
x = symbols("x")
init_printing()
print(latex(Integral(sqrt(1/x), x)))
```

`## \int \sqrt{\frac{1}{x}}\, dx`