Direct & Iterative Methods for Solving Linear Systems: Gauss, Crout, SOR, Gauss-Seidel

TAKEAWAYS:

• Direct methods (Gauss, Crout) give exact solutions in finite steps but are computationally expensive for large systems.
• Iterative methods (Jacobi, Gauss-Seidel, SOR) approximate solutions and are efficient for sparse or large systems, but require convergence analysis.
• The Gauss-Seidel method converges faster than Jacobi because it uses updated values immediately.
• Successive Over-Relaxation (SOR) improves convergence by introducing a relaxation factor ω (1 < ω < 2).
• Error bounds and residuals determine method accuracy; iterative methods stop when residuals fall below a tolerance.
• Real-world applications include electric circuit analysis (Nepal Electricity Authority’s load distribution), financial modeling (bank loan repayment schedules), and traffic routing (Pathao’s dynamic ride-matching).

1. Introduction: Why Solve Linear Systems Numerically?

Linear systems arise in engineering when modeling equilibrium states, such as:

• Structural analysis: Forces in a bridge truss (Nepal’s Seti River Bridge).
• Circuit analysis: Voltage/current in NTC’s power grids.
• Economics: Input-output models in Nepal’s GDP calculations.

Key Definitions:

• Linear system: $A\mathbf{x}=\mathbf{b}$, where $A$ is an $n\times n$ matrix, $\mathbf{x}$ is the solution vector, and $\mathbf{b}$ is the right-hand side.
• Direct methods: Compute exact solutions in finite steps (e.g., Gauss elimination).
• Iterative methods: Approximate solutions via successive refinements (e.g., Jacobi, Gauss-Seidel).

2. Direct Methods: Exact Solutions

2.1 Gauss Elimination Method

How it works:

1. Convert the system into row-echelon form (upper triangular) using row operations.
2. Perform back substitution to solve for variables.

Example 1: Solve using Gauss elimination: $$\{\begin{array}{l}3x-y+z=10\\ 2x+3y+2z=14\\ x+2y+3z=14\end{array}$$

Step-by-Step Trace:

1. Eliminate $x$ from rows 2 and 3:
   • Row2 = Row2 – (2/3)Row1
   • Row3 = Row3 – (1/3)Row1
2. Eliminate $y$ from row 3:
   • Row3 = Row3 – (7/11)Row2
3. Back substitution:
   • $z=0$
   • $y=(8/3-(4/3)(0))/(11/3)=8/11$
   • $x=(10-(-1)(8/11)+0)/3=100/33$

Solution: $x=\frac{100}{33}$, $y=\frac{8}{11}$, $z=0$.

Advantages:

• Exact solution (no rounding errors in theory).
• Simple to implement.

Disadvantages:

• Computationally expensive for large $n$ ($O(n^3)$ operations).
• Prone to round-off errors in floating-point arithmetic.

2.2 Crout Factorization (LU Decomposition)

How it works: Decompose $A=LU$, where:

• $L$: Lower triangular matrix with 1s on the diagonal.
• $U$: Upper triangular matrix.

Example 2: Solve using Crout’s method: $$\{\begin{array}{l}x+4y+7z=4\\ 4x+4y+3z=8\\ x+6y+2z=6\end{array}$$

Step-by-Step Trace:

1. Factorize $A$: $$A=\left[\begin{array}{ccc}1& 4& 7\\ 4& 4& 3\\ 1& 6& 2\end{array}\right]=LU$$ Solve for $L$ and $U$: $$L=\left[\begin{array}{ccc}1& 0& 0\\ 4& 1& 0\\ 1& l_{32}& 1\end{array}\right],U=\left[\begin{array}{ccc}1& 4& 7\\ 0& u_{22}& u_{23}\\ 0& 0& u_{33}\end{array}\right]$$ Using $A=LU$, compute:

   • $u_{22}=4-4\times 4=-12$
   • $l_{32}=(6-4\times 1)/(-12)=-1/3$
   • $u_{33}=2-7\times 1-3\times (-1/3)=-2/3$

2. Solve $L\mathbf{y}=\mathbf{b}$: $$\left[\begin{array}{ccc}1& 0& 0\\ 4& 1& 0\\ 1& -1/3& 1\end{array}\right]\left[\begin{array}{c}{y}_1\\ y_2\\ y_3\end{array}\right]=\left[\begin{array}{c}4\\ 8\\ 6\end{array}\right]$$ Solution: $y_1=4$, $y_2=-8$, $y_3=10$.

3. Solve $U\mathbf{x}=\mathbf{y}$: $$\left[\begin{array}{ccc}1& 4& 7\\ 0& -12& 3\\ 0& 0& -2/3\end{array}\right]\left[\begin{array}{c}{x}_1\\ x_2\\ x_3\end{array}\right]=\left[\begin{array}{c}4\\ -8\\ 10\end{array}\right]$$ Solution: $x_3=-15$, $x_2=1$, $x_1=-1$.

Solution: $x=-1$, $y=1$, $z=-15$.

Advantages:

• Efficient for repeated solutions with the same $A$ (e.g., optimization problems).
• Avoids pivoting issues if $A$ is well-conditioned.

Disadvantages:

• Requires $O(n^3)$ operations.
• Fails if $A$ is singular or ill-conditioned.

3. Iterative Methods: Approximate Solutions

3.1 Jacobi Method

How it works: Rewrite $A\mathbf{x}=\mathbf{b}$ as: $${\mathbf{x}}^{(k+1)}=D^{-1}(-E{\mathbf{x}}^{(k)}-F{\mathbf{x}}^{(k)}+\mathbf{b})$$ where $A=D+E+F$, $D$ is diagonal, $E$ is lower triangular, and $F$ is upper triangular.

Convergence Criterion: Converges if $\rho (D^{-1}(E+F))<1$, where $\rho$ is the spectral radius.

Example 3: Apply Jacobi to: $$\{\begin{array}{l}4x-y+z=4\\ x+5y+z=7\\ x+y+3z=3\end{array}$$ Initial guess: ${\mathbf{x}}^{(0)}=[0,0,0{]}^T$.

Iteration Steps:

flowchart TD
    A["Iteration 0: x₀=0, y₀=0, z₀=0"] --> B["Compute x₁ = (4 + y₀ - z₀)/4 = 1"]
    B --> C["Compute y₁ = (7 - x₁ - z₀)/5 = 6/5 = 1.2"]
    C --> D["Compute z₁ = (3 - x₁ - y₁)/3 ≈ 0.133"]
    D --> E["Iteration 1: x₁=1, y₁=1.2, z₁≈0.133"]
    E --> F["Repeat until convergence (|x^(k+1) - x^(k)| < ε)"]

After 5 iterations: $\mathbf{x}\approx [0.999,1.001,0.999{]}^T$.

Advantages:

• Simple to implement.
• Parallelizable (each component updated independently).

Disadvantages:

• Slow convergence (often requires many iterations).
• May not converge for all systems.

3.2 Gauss-Seidel Method

Improvement over Jacobi: Uses the most recent values immediately.

Example 4: Solve the same system as Example 3 using Gauss-Seidel.

Iteration Steps:

flowchart TD
    A["Iteration 0: x₀=0, y₀=0, z₀=0"] --> B["Compute x₁ = (4 + y₀ - z₀)/4 = 1"]
    B --> C["Compute y₁ = (7 - x₁ - z₀)/5 = 6/5 = 1.2"]
    C --> D["Compute z₁ = (3 - x₁ - y₁)/3 ≈ 0.133"]
    D --> E["Iteration 1: x₁=1, y₁=1.2, z₁≈0.133"]
    E --> F["Use updated x₁, y₁ for next iteration"]
    F --> G["Iteration 2: x₂ ≈ 0.999, y₂ ≈ 1.000, z₂ ≈ 0.999"]

Converges in fewer iterations than Jacobi.

Convergence Criterion: Converges if $A$ is strictly diagonally dominant or symmetric positive definite.

3.3 Successive Over-Relaxation (SOR)

How it works: Introduce a relaxation factor $\omega$ (typically $1<\omega <2$) to accelerate convergence: $${\mathbf{x}}^{(k+1)}=(1-\omega ){\mathbf{x}}^{(k)}+\omega D^{-1}(-E{\mathbf{x}}^{(k+1)}-F{\mathbf{x}}^{(k)}+\mathbf{b})$$

Example 5: Solve Example 3 with $\omega =1.25$.

Optimal $\omega$: For symmetric positive definite matrices, ${\omega }_{\text{opt}}=\frac{2}{1+\sqrt{1-{\rho }^2}}$, where $\rho$ is the spectral radius of $D^{-1}(E+F)$.

Comparison Table:

4. Error Analysis and Stopping Criteria

Error Measures:

1. Residual: ${\mathbf{r}}^{(k)}=\mathbf{b}-A{\mathbf{x}}^{(k)}$. Stop when $\Vert {\mathbf{r}}^{(k)}\Vert <ϵ$.
2. Relative Error: $\frac{\Vert {\mathbf{x}}^{(k+1)}-{\mathbf{x}}^{(k)}\Vert }{\Vert {\mathbf{x}}^{(k+1)}\Vert }<ϵ$.

Example 6: For Example 3, stop when $\Vert {\mathbf{r}}^{(k)}\Vert <10^{-3}$.

5. Real-World Applications

5.1 Nepal Electricity Authority (NEA) Load Distribution

• Problem: Distribute power across Nepal’s grid to minimize losses.
• Method: Gauss-Seidel solves the linear system $A\mathbf{V}=\mathbf{I}$, where $\mathbf{V}$ is voltage and $\mathbf{I}$ is current.
• Why? Sparse matrix $A$ (only nearby nodes interact), so iterative methods are efficient.

5.2 Bank Loan Repayment (Nabil Bank)

• Problem: Calculate monthly payments for a loan with interest.
• Method: Solve $A\mathbf{x}=\mathbf{b}$, where $A$ encodes interest rates and $\mathbf{x}$ is the repayment schedule.
• Why? Crout factorization is used for repeated calculations (e.g., adjusting interest rates).

5.3 Pathao’s Dynamic Ride-Matching

• Problem: Assign drivers to passengers in real-time to minimize wait times.
• Method: Gauss-Seidel approximates optimal assignments by iteratively improving matches.
• Why? Large-scale system (thousands of variables), so iterative methods outperform direct methods.

6. Exam Tip: What to Expect

1. Direct Methods (Gauss, Crout):

   • Expect step-by-step elimination or factorization questions.
   • Common pitfall: Forgetting to perform back substitution after elimination.
   • Tip: Always write the augmented matrix clearly.

2. Iterative Methods (Jacobi, Gauss-Seidel, SOR):

   • Always check convergence before submitting answers.
   • For SOR: State the value of $\omega$ used (e.g., $\omega =1.25$).
   • Tip: Show the first 2–3 iterations explicitly.

3. Error Analysis:

   • Questions may ask for residuals or relative error after a fixed number of iterations.
   • Tip: Use $\Vert \mathbf{r}\Vert =\Vert \mathbf{b}-A\mathbf{x}\Vert$ for residuals.

4. Matrix Properties:

   • Diagonally dominant matrices guarantee Gauss-Seidel convergence.
   • Symmetric positive definite matrices have optimal $\omega$ for SOR.

5. Past Exam Patterns:

   • Part (a): Direct method (Gauss/Crout) with 3 equations.
   • Part (b): Iterative method (Jacobi/Gauss-Seidel/SOR) with 3–4 equations.
   • Always verify your solution by plugging back into the original equations.

7. Summary of Key Formulas

8. Worked Example: Traffic Flow Optimization (Kathmandu)

Problem: Model traffic flow at a 3-intersection junction in Kathmandu using the system: $$\{\begin{array}{l}2x-y+z=100\\ -x+3y+z=150\\ x+y+2z=200\end{array}$$ where $x,y,z$ are traffic flows (vehicles/hour).

Solution:

1. Gauss Elimination:

   • Eliminate $x$ from rows 2 and 3: $$\left[\begin{array}{ccccc}2& -1& 1& |& 100\\ 0& 2.5& 1.5& |& 200\\ 0& 0.5& 1.5& |& 150\end{array}\right]$$
   • Back substitution yields: $z=50$, $y=50$, $x=50$.

2. Gauss-Seidel Iteration:

   • Initial guess: ${\mathbf{x}}^{(0)}=[0,0,0{]}^T$.
   • After 3 iterations: ${\mathbf{x}}^{(3)}\approx [50.0,50.0,50.0{]}^T$.

Interpretation:

• Equal flow ($x=y=z=50$ vehicles/hour) optimizes traffic at the junction.
• Real-world tie-in: Used by SmartRide (Nepal’s traffic management app) to simulate and suggest reroutes.

9. Common Mistakes to Avoid

1. Floating-Point Errors:

   • In Gauss elimination, partial pivoting (swap rows to avoid division by small numbers) is often required.
   • Example: Swap rows if $|a_{ij}|<ϵ$.

2. Incorrect Iterative Updates:

   • In Gauss-Seidel, always use the latest values for $j<i$.
   • Mistake: Using old values for all $j\ne i$ (this is Jacobi, not Gauss-Seidel).

3. Forgetting Convergence Checks:

   • Always state the stopping criterion (e.g., $\Vert \mathbf{r}\Vert <10^{-4}$).

4. Misapplying SOR:

   • $\omega$ must be chosen carefully. For Example 3, $\omega =1.25$ works, but $\omega =1.5$ may diverge.

10. Practice Questions (Exam-Style)

1. Gauss Elimination: Solve: $$\{\begin{array}{l}2x+y-z=8\\ -3x-y+2z=-11\\ -2x+y+2z=-3\end{array}$$

2. Crout Factorization: Factorize and solve: $$\{\begin{array}{l}4x+y+z=4\\ 3x+4y+z=5\\ x+y+4z=6\end{array}$$

3. Gauss-Seidel: Apply Gauss-Seidel to the above system with ${\mathbf{x}}^{(0)}=[0,0,0{]}^T$. Stop when $\Vert \mathbf{r}\Vert <0.01$.

4. SOR: Solve the same system with $\omega =1.1$. Compare convergence speed to Gauss-Seidel.

11. Further Reading

• Books:
  • Numerical Methods for Engineers by Steven C. Chapra.
  • Introduction to Numerical Analysis by Kendall E. Atkinson.
• Online Resources:
  • MIT OpenCourseWare: Numerical Linear Algebra.
  • Python libraries: numpy.linalg.solve (direct), scipy.sparse.linalg (iterative).