Jacobi Method in Python and Numpy

ln this, we will discuss the Jacobi Method implementation in Python. The Jacobi method (or the Jacobi iterative method is an algorithm for determining solutions for a system of linear equations that diagonally dominant.

Jocobi Method

Jacobi method is a matrix iterative method used to solve the linear equation Ax = b of a known square matrix of magnitude n * n and vector b or length n.

Jacobi's method is widely used in boundary calculations (FDM), which is an important part of the financial world. It can be done in such a way that it is solved by finite difference technique. Jacobi method is one of the ways to solve the resulting matrix equation that arises from FDM.

Jocobi Method with Numpy

We will use the NumPy library to speed up the calculation of the Jacobi method. NumPy works much better than writing implementations in pure Python. 

The iterative nature of the Jacobi method means that any increase in speed within each iteration can have a significant impact on the overall calculation.

Now take a look at how to implement Jocobi Method in Python Programming using Numpy library.

Example 1: Implement of the following 

A = [[ 4. -2.  1.]    ,   b [1.0, 2.0, 3.0]    , x [1.0, 1.0, 1.0]
       [ 1. -3.  2.]
       [-1.  2.  6.]]

Python code:

import numpy as np
from scipy.linalg import solve
def jacobi(A,b,x,n):
    D = np.diag(A)
    R = A-np.diagflat(D)
    for i in range(n):
        x = (b-np.dot(R,x))/D
    return x
A = np.array([[4,-2.0,1.0],[1.0,-3.0,2.0],[-1.0,2.0,6.0]])
b = [1.0,2.0,3.0]
x = [1.0,1.0,1.0]
n = 25
x = jacobi(A,b,x,n)
print(solve(A,b))

Output:

A [[ 4. -2.  1.]
 [ 1. -3.  2.]
 [-1.  2.  6.]]
b [1.0, 2.0, 3.0]
x [1.0, 1.0, 1.0]

[-0.04109589 -0.28767123  0.5890411 ]

Example 2: Implement of the following

A= [[2 1] , b= [11, 13] , x= [1.0, 1.0]
[5 7]]

Python Code:

import numpy as np
from scipy.linalg import solve
def jacobi(A,b,x,n):
    D = np.diag(A)
    R = A-np.diagflat(D)
    for i in range(n):
        x = (b-np.dot(R,x))/D
    return x
A = np.array([[2,1],[5,7]])
print(';A=',A)
b = [11,13]
print('b=',b)
x = [1.0,1.0]
print('x=',x)
n = 25
x = jacobi(A,b,x,n)
print(solve(A,b))

Output:

A= [[2 1]
 [5 7]]
b= [11, 13]
x= [1.0, 1.0]
[ 7.11111111 -3.22222222]