Solving a system of linear equations using Jacobi's method

jacobi.py

# -*- coding: utf-8 -*-
"""
Created on Sun Jan 19 11:53:03 2020

@author: DASA0
"""

import numpy as np

a = np.array([[1, 2, -3, 4],
           [2, 2, -2, 3],
           [0, 1, 1, 0],
           [1, -1, 1, -2]])
b = np.array([12, 10, -1, -4])
(n,) = np.shape(b)
x = np.full(n, 1.0, float) # initial guess
xnew = np.empty(n, float)
iterlimit = 100
tolerance = 1.0e-4

# iterations:
for iter in range(iterlimit):
    for i in range(n):
        s = 0
        for j in range(n):
            if j != i:
                s += a[i, j]*x[j]
        xnew[i] = -1/a[i,i] * (s - b[i])
    if (abs(xnew - x) < tolerance).all():
        break
    else:
        x = np.copy(xnew)

print('Number of iterations: %d '% (iter + 1))
print('The solution for the system of linear equations is:')
print(x)

The above example also shows that since the diagonal dominance condition is not satisfied, the results obtained are vague. For such cases , Gauss elimination is preferred to iterative methods such as Jacobi or Gauss-Siedel.