SoFunction
Updated on 2024-11-14

Detailed Python Implementation of Jacobi Iteration Algorithm

import numpy as np
import time

1.1 Jacobi Iterative Algorithm

def Jacobi_tensor_V2(A,b,Delta,m,n,M):
start=time.perf_counter()#Start the clock
find=0# Used to mark if convergence occurs within a specified number of steps
X=(n)#Iteration start point
x=(n)# Used to store intermediate results of iterations
d=(n)# for storing the diagonal portion of Ax**(m-2)
m1=m-1
m2=2-m
for i in range(M):
print('X',X)
a=(A)
# Got Ax**(m-2)
for j in range(m-2):
a=(a,X)
#Get d and (2-m)Dx**(m-2)+(L'+U')x**(m-2)
for j in range(n):
d[j]=a[j,j]
a[j,j]=m2*a[j,j]
#Iterative updates
for j in range(n):
x[j]=(b[j]-(a[j],X))/(m1*d[j])
#Determine if accuracy requirements are met
if ((X-x))<Delta:
find=1
break 
X=(x)
end=time.perf_counter()# End of the clock
print('Time:',end-start)
print('Iteration',i)
return X,find,i,end-start

1.2 Generating function for tensor A and generating function for vector b:

def Creat_A(m,n):# Generate a tensor A
size=(m, n)
X=(n)
while 1:
# Randomly generate a tensor A of a given shape
A=(-49,50,size=size)
# Determine if Dx**(m-2) is non-singular, if so, the requirement is satisfied and jump out of the loop
D=(A)
for i1 in range(n):
for i2 in range(n):
if i1!=i2:
D[i1,i2]=0
for i in range(m-2):
D=(D,X)
det=(D)
if det!=0:
break
# Expand the diagonal tensor of A by a factor of ten so that the diagonal is dominant
for i1 in range(n):
for i2 in range(n):
if i1==i2:
A[i1,i2]=A[i1,i2]*10
print('A:')
print(A)
return A
# Generate vector b from A and given X according to Ax**(m-1)=b
def Creat_b(A,X,m):
a=(A)
for i in range(m-1):
a=(a,X)
print('b:')
print(a)
return a

1.3 Generating function for the symmetric tensor S:

def Creat_S(m,n):# Generate a symmetric tensor B
size=(m, n)
S=(size)
print('S',S)
for i in range(4):
# Generate n as a vector a
a=(n)*(-5,6)
b=(a)
# Perform m-1 outer products on a to obtain a rank 1 symmetric tensor b
for j in range(m-1):
b=outer(b,a)
# Stacking different b's yields the low-rank symmetric tensor S
S=S+b
print('S:')
print(S)
return S
def outer(a,b):
c=[]
for i in b:
(i*a)
return (c)
return a

1.4 Experiment I

def test_1():
Delta=0.01# Accuracy
m=3# Order of A
n=3# Dimension of A
M=200# Maximum number of iteration steps
X_real=( [2,3,4])
A=Creat_A(m,n) 
b=Creat_b(A,X_real,m)
Jacobi_tensor_V2(A,b,Delta,m,n)

This is the whole content of this article.