Implementation of Jacobi Method to Solve a System of Linear Equations in Python

The Jacobi Method is an iterative algorithm for solving systems of linear equations. It starts with initial guesses and repeatedly refines them until convergence is achieved.

Mathematical Foundation

Consider a system of linear equations:

a??x? + a??x? + ... + a??x? = b?

a??x? + a??x? + ... + a??x? = b?

...

a??x? + a??x? + ... + a??x? = b?

The Jacobi method rearranges each equation to isolate one variable. For the i-th equation:

x_i = (b_i - ?(a_ijx_j)) / a_ii

The Jacobi Algorithm

The algorithm follows these steps:

• Start with initial guesses for all unknowns

• Calculate new values using the rearranged equations

• Check if the error between new and old values is below the tolerance

• If not converged, update guesses and repeat

Example Problem

Let's solve this system:

20x + y - 2z = 17

3x + 20y - z = -18

2x - 3y + 20z = 25

Rearranging for each variable:

x = (17 - y + 2z) / 20

y = (-18 - 3x + z) / 20

z = (25 - 2x + 3y) / 20

Implementation Using Individual Equations

import numpy as np

# Initial guesses
x, y, z = 0.0, 0.0, 0.0
error = 1.0
tolerance = 1e-5
iteration = 0
max_iterations = 100

print(f"{'Iteration':<10}{'x':<12}{'y':<12}{'z':<12}{'Error':<12}")
print("-" * 60)

while error > tolerance and iteration < max_iterations:
    iteration += 1
    
    # Store old values
    x_old, y_old, z_old = x, y, z
    
    # Calculate new values
    x = (17 - y_old + 2*z_old) / 20
    y = (-18 - 3*x_old + z_old) / 20
    z = (25 - 2*x_old + 3*y_old) / 20
    
    # Calculate error
    error = abs(x - x_old) + abs(y - y_old) + abs(z - z_old)
    
    print(f"{iteration:<10}{x:<12.6f}{y:<12.6f}{z:<12.6f}{error:<12.6f}")

print(f"\nSolution: x = {x:.6f}, y = {y:.6f}, z = {z:.6f}")

Iteration x           y           z           Error       
------------------------------------------------------------
1         0.850000    -0.900000   1.250000    3.000000    
2         1.035000    -1.022500   0.992500    0.357500    
3         0.996125    -0.998875   1.001625    0.062500    
4         1.000406    -0.999949   0.999719    0.006219    
5         0.999968    -1.000008   1.000016    0.000476    
6         1.000003    -1.000000   0.999998    0.000036    
7         1.000000    -1.000000   1.000000    0.000003    

Solution: x = 1.000000, y = -1.000000, z = 1.000000

Matrix Implementation

A more general approach uses matrix operations ?

import numpy as np

def jacobi_method(A, b, x0=None, tolerance=1e-5, max_iterations=100):
    """
    Solve Ax = b using Jacobi method
    """
    n = len(b)
    if x0 is None:
        x = np.zeros(n)
    else:
        x = x0.copy()
    
    x_new = np.zeros(n)
    
    print(f"{'Iteration':<10}{'Solution Vector':<30}{'Error':<12}")
    print("-" * 55)
    
    for iteration in range(max_iterations):
        # Calculate new values
        for i in range(n):
            sum_ax = sum(A[i][j] * x[j] for j in range(n) if i != j)
            x_new[i] = (b[i] - sum_ax) / A[i][i]
        
        # Calculate error
        error = np.sum(np.abs(x_new - x))
        
        print(f"{iteration+1:<10}{str(np.round(x_new, 6)):<30}{error:<12.6f}")
        
        if error < tolerance:
            break
            
        x = x_new.copy()
    
    return x_new, iteration + 1

# Define 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)

# Solve the system
solution, iterations = jacobi_method(A, b)

print(f"\nFinal solution: {solution}")
print(f"Converged in {iterations} iterations")

Iteration Solution Vector               Error       
-------------------------------------------------------
1         [ 0.85 -0.9   1.25]           3.000000    
2         [ 1.035    -1.0225   0.9925 ] 0.357500    
3         [ 0.996125 -0.998875  1.001625] 0.062500    
4         [ 1.000406 -0.999949  0.999719] 0.006219    
5         [ 0.999968 -1.000008  1.000016] 0.000476    
6         [ 1.000003 -1.       0.999998] 0.000036    
7         [ 1. -1.  1.]                 0.000003    

Final solution: [ 1. -1.  1.]
Converged in 7 iterations

Key Points

• The method requires the diagonal elements to be non-zero and dominant

• Convergence depends on the matrix properties diagonally dominant matrices converge faster

• Each iteration uses values from the previous iteration only

• The method is suitable for sparse matrices and parallel computation

Comparison with Other Methods

Conclusion

The Jacobi method is an effective iterative technique for solving linear systems, especially for large sparse matrices. While it converges slower than Gauss-Seidel, its parallelizable nature makes it valuable for modern computing applications.