Modelling the Gauss Seidel Method in Python

The Gauss-Seidel method is an iterative technique for solving systems of linear equations. Unlike the Jacobi method, Gauss-Seidel uses newly computed values within the same iteration, which speeds up convergence.

Mathematical Foundation

A system of linear equations can be written as:

$$\mathrm{a_{1,1}x_{1} \: + \: a_{1,2}x_{2} \: + \: \dotso \: + \: a_{1,n}x_{n} \: = \: b_{1}}$$

$$\mathrm{a_{2,1}x_{1} \: + \: a_{2,2}x_{2} \: + \: \dotso \: + \: a_{2,n}x_{n} \: = \: b_{2}}$$

$$\mathrm{\vdots}$$

$$\mathrm{a_{n,1}x_{1} \: + \: a_{n,2}x_{2} \: + \: \dotso \: + \: a_{n,n}x_{n} \: = \: b_{n}}$$

The general Gauss-Seidel formula is:

$$\mathrm{x_{i_{new}} \: = \: - \frac{1}{a_{i,i}}(\sum_{j=1,j \neq i}^n a_{i,j}x_{j_{new}} \: - \: b_{i})}$$

where $\mathrm{x_{j_{new}}}$ represents newly computed values from the current iteration.

Algorithm Steps

• Start with initial guesses for all unknowns

• For each equation, calculate new values using the rearranged form

• Use newly computed values immediately in subsequent calculations

• Calculate error: $\mathrm{|x_{new} \: - \: x_{guess}|}$

• If error > tolerance (e.g., $\mathrm{10^{-5}}$), continue iterating

• Otherwise, the solution has converged

Convergence Criterion

For guaranteed convergence, the system must satisfy diagonal dominance (Scarborough criterion):

$$\mathrm{\frac{\sum_{i\neq j |a_{i,j}|}}{|a_{i,i}|} \begin{cases} \le \: 1 \: \text{for all equations} \

Example Problem

Let's solve this system of equations:

$$\mathrm{20x \: + \: y \: - \: 2z \: = \: 17}$$

$$\mathrm{3x \: + \: 20y \: - \: z \: = \: -18}$$

$$\mathrm{2x \: - \: 3y \: + \: 20z \: = \: 25}$$

Rearranging for each variable:

$$\mathrm{x_{new} \: = \: (-y_{guess} \: + \: 2z_{guess} \: + \: 17)/20}$$

$$\mathrm{y_{new} \: = \: (-3x_{new} \: + \: z_{guess} \: - \: 18)/20}$$

$$\mathrm{z_{new} \: = \: (-2x_{new} \: + \: 3y_{new} \: + \: 25)/20}$$

Method 1: Equation-by-Equation Implementation

import numpy as np
import matplotlib.pyplot as plt

# Initial guess
x, y, z = 0, 0, 0
error = 1
count = 0
errors = []

while error > 1e-5:
    count += 1
    
    # Store previous values
    x_prev, y_prev, z_prev = x, y, z
    
    # Update using newly computed values
    x = (17 - y + 2*z) / 20
    y = (z - 18 - 3*x) / 20     # Uses new x
    z = (25 - 2*x + 3*y) / 20   # Uses new x and y
    
    # Calculate error
    error = abs(x - x_prev) + abs(y - y_prev) + abs(z - z_prev)
    errors.append(error)

print(f"Solution: x = {x:.5f}, y = {y:.5f}, z = {z:.5f}")
print(f"Converged in {count} iterations")

Solution: x = 1.00000, y = -1.00000, z = 1.00000
Converged in 5 iterations

Method 2: Matrix Implementation

For larger systems, a matrix approach is more practical:

import numpy as np

# Coefficient matrix and RHS vector
A = np.array([[20, 1, -2],
              [3, 20, -1],
              [2, -3, 20]], dtype=float)
b = np.array([17, -18, 25], dtype=float)

n = len(b)
x = np.zeros(n)  # Initial guess
error = 1
count = 0

while error > 1e-5:
    count += 1
    x_prev = x.copy()
    
    for i in range(n):
        sum_ax = 0
        for j in range(n):
            if i != j:
                sum_ax += A[i, j] * x[j]
        
        x[i] = (b[i] - sum_ax) / A[i, i]
    
    # Calculate error
    error = np.sum(np.abs(x - x_prev))

print(f"Solution: {x}")
print(f"Converged in {count} iterations")

Solution: [ 1. -1.  1.]
Converged in 5 iterations

Key Advantages

• Faster convergence: Uses updated values immediately

• Memory efficient: Only stores one set of variables

• Simple implementation: Straightforward algorithm

Conclusion

The Gauss-Seidel method efficiently solves systems of linear equations through iterative refinement. Success depends on diagonal dominance ensure the coefficient matrix satisfies this condition for guaranteed convergence.