Python 线性方程组求解之：Jacobi迭代算法

1. import numpy as np
2. delta=0.0001#精度要求
3. #数据读取
4. data = []
5. f = open("h:/notepad/数学实验/jacobi_data.txt")
6. for line in f:
7. line = line.replace("\n","")
8. data.append(list(map(eval, line.split(","))))
9. f.close()
10. data=np.array(data)
11. #对data行初等变换，主元变为每列最大值
12. row,column=data.shape
13. for i in range(row):
14. max_value_index=np.argmax(np.fabs(data[i:row,i]))
15. temp=np.copy(data[i,:])
16. data[i,:]=data[max_value_index+i,:]
17. data[max_value_index+i,:]=temp
18. #lu: -(l+u)
19. #d:系数矩阵的对角线元素
20. #b:ax=b中的b
21. lu=np.negative(data[:,0:column-1])
22. d=np.zeros(row)
23. b=data[:,column-1]
24. for i in range(row):
25. d[i]=data[i,i]
26. lu[i,i]=0
28. #迭代求解
29. x=np.ones(row)#用于存储迭代过程中x的值
30. y=np.ones(row)#用于存储中间结果
31. for iteration in range(100):
32. print('x:',x)
33. #迭代计算
34. for i in range(row):
35. y[i]=np.vdot(lu[i,:],x)+b[i]
36. y[i]=y[i]/d[i]
38. #判断是否达到精度要求
39. if np.max(np.fabs(x-y))<delta:
40. print('iteration:',iteration)
41. break
42. #将y幅值到x，开始下一轮迭代
43. x=np.copy(y)

原理：

注：

 

实例数据：jacobi_data.txt 。组织方式：[a,b]

2,-1,-1,-5
1,1,10,11
1,5,-1,8