Simultaneous Linear Equations

Created: 2019-12-11 ᐱ 15:32

Learning Objectives

  • Use the Gaussian Elimination method of solving sets of simultaneous linear equations.
  • Select the correct ordering of rows to minimize floating point error.

The Problem

  • We are given a series of equations of this form:

\[ \begin{array}{l} a_1x + b_1y + c_1z = d_1 \\ a_2x + b_2y + c_2z = d_2 \\ a_3x + b_3y + c_3z = d_3 \\ \end{array} \]

  • We know the values of \(a_i,~b_i,~\hbox{and}~c_i\) for all \(i\)
  • We want the values of \(x,y,~\hbox{and}~z\).
  • Solution: use a matrix formulation
\[ \begin{array}{llll} \left[\begin{array}{lll} a_1 & b_1 & c_1 \\ a_2 & b_2 & c_2 \\ a_3 & b_3 & c_3 \\ \end{array}\right] & \left[\begin{array}{l} x \\ y \\ z \\ \end{array}\right] & = & \left[\begin{array}{l} d_1 \\ d_2 \\ d_3 \\ \end{array}\right] \end{array} \]

Solving The matrix equations

  • Use Gaussian Elimination to solve this. Step 1: Forward elimination.
    • For each row \(i\) and row \(j\) where \(i \lt j\), subtract off multiples of row \(i\) from row \(j\) to zero out column \(i\).
    • The solution vector needs to participate in this as well: called the augmented matrix
\[ \begin{array}{lllll} \left[\begin{array}{lll|l} a_1 & b_1 & c_1 & d_1 \\ a_2 & b_2 & c_2 & d_2 \\ a_3 & b_3 & c_3 & d_3 \\ \end{array}\right] & \Rightarrow & \left[\begin{array}{lll|l} a_1 & b_1 & c_1 & d_1 \\ 0 & b'_2 & c'_2 & d'_2 \\ 0 & b'_3 & c'_3 & d'_3 \\ \end{array}\right] & \Rightarrow & \left[\begin{array}{lll|l} a_1 & b_1 & c_1 & d_1 \\ 0 & b'_2 & c'_2 & d'_2 \\ 0 & 0 & c''_3 & d''_3 \\ \end{array}\right] \end{array} \]
  • Step 2: Backward elimination.
    • For each row \(i\) and row \(j\) where \(i>j\), subtract off multiples of row \(i\) from row \(j\) to zero out column \(j\).
    • Result: a diagonal matrix. Then \(z = d''_3 / c''_3\), etc…

A Real Example, Step 1

  • Want to solve:
\[ \begin{array}{llll} \left[\begin{array}{lll} 1 & 3 & 5 \\ 2 & 10 & 30 \\ 3 & 12 & 20 \\ \end{array}\right] & \left[\begin{array}{l} x \\ y \\ z \\ \end{array}\right] & = & \left[\begin{array}{l} 10 \\ 30 \\ 50 \\ \end{array}\right] \end{array} \]
  • Step 1: forward substitution
\[ \begin{array}{lllll} \left[\begin{array}{lll|l} 1 & 3 & 5 & 10 \\ 2 & 10 & 30 & 30 \\ 3 & 17 & 20 & 50 \\ \end{array}\right] & \Rightarrow & \left[\begin{array}{lll|l} 1 & 3 & 5 & 10 \\ 0 & 4 & 20 & 10 \\ 0 & 8 & 5 & 20 \\ \end{array}\right] & \Rightarrow & \left[\begin{array}{lll|l} 1 & 3 & 5 & 10 \\ 0 & 4 & 20 & 10 \\ 0 & 0 & -35 & 0 \\ \end{array}\right] \end{array} \]

Step 2

\[ \begin{array}{llll} \left[\begin{array}{lll} 1 & 3 & 5 \\ 2 & 10 & 30 \\ 3 & 12 & 20 \\ \end{array}\right] & \left[\begin{array}{l} x \\ y \\ z \\ \end{array}\right] & = & \left[\begin{array}{l} 10 \\ 30 \\ 50 \\ \end{array}\right] \end{array} \]
  • Step 2: Backward substitution
\[ \begin{array}{lllll} \left[\begin{array}{lll|l} 1 & 3 & 5 & 10 \\ 0 & 4 & 20 & 10 \\ 0 & 0 & -35 & 0 \\ \end{array}\right] & \Rightarrow & \left[\begin{array}{lll|l} 1 & 3 & 0 & 10 \\ 0 & 4 & 0 & 10 \\ 0 & 0 & -35 & 0 \\ \end{array}\right] & \Rightarrow & \left[\begin{array}{lll|l} 1 & 0 & 0 & 2.5 \\ 0 & 4 & 0 & 10 \\ 0 & 0 & -35 & 0 \\ \end{array}\right] \end{array}\]
  • Solution: \( x=2.5, y=1.25, z = 0 \)
  • Tips:
    • If you end up with a all-zero row, the system is underspecified.
    • To reduce numerical error, sort rows by largest column value.

Source Code

 1: #define MAX_N 100
 2: struct AugmentedMatrix { double mat[MAX_N][MAX_N + 1]; };
 3: struct ColumnVector { double vec[MAX_N]; };
 4: ColumnVector GaussianElimination(int N, AugmentedMatrix Aug) {// O(N^3)
 5:     // input: N, Augmented Matrix Aug, output: Column vector X, the answer
 6:     int i, j, k, l; double t; ColumnVector X;
 7:     // the forward elimination phase
 8:     for (j = 0; j < N - 1; j++) {
 9:        l = j;
10:        for (i = j + 1; i < N; i++) // which row has largest column value
11:           if (fabs(Aug.mat[i][j]) > fabs(Aug.mat[l][j])) l = i;
12:        // swap this pivot row, reason: to minimize floating point error
13:        for (k = j; k <= N; k++)
14:            t = Aug.mat[j][k], Aug.mat[j][k] = Aug.mat[l][k], Aug.mat[l][k] = t;
15:         // the actual forward elimination phase
16:        for (i = j + 1; i < N; i++)
17:            for (k = N; k >= j; k--)
18:                Aug.mat[i][k] -= Aug.mat[j][k] * Aug.mat[i][j] / Aug.mat[j][j];
19:     }
20:     for (j = N - 1; j >= 0; j--) {// the back substitution phase
21:        for (t = 0.0, k = j + 1; k < N; k++) t += Aug.mat[j][k] * X.vec[k];
22:        X.vec[j] = (Aug.mat[j][N] - t) / Aug.mat[j][j];
23:     }
24:     return X;
25: }