用Python實現雅可比方法求解線性方程組

---

這是解決如下所示線性方程組最直接的迭代方法。

$$\mathrm{a_{1,1}\: x_{1} \: + \: a_{1,2} \: x_{2} \: + \: \dotso\dotso \: + \: a_{1,n} \: x_{n} \: = \: b_{1}}$$

$$\mathrm{a_{2,1} \: x_{1} \: + \: a_{2,2} \: x_{2} \: + \: \dotso\dotso \: + \: a_{2,n} \: x_{n} \: = \: b_{2}}$$

$$\mathrm{\vdots}$$

$$\mathrm{a_{n,1} \: x_{1} \: + \: a_{n,2} \: x_{2} \: + \: \dotso\dotso \: + \: a_{n,n} \: x_{n} \: = \: b_{n}}$$

基本思想是：將每個線性方程重新排列，將一個新的變數移到左邊。然後，從對每個變數的初始猜測開始確定新值。在接下來的迭代中，新值將用作更準確的猜測。重複此迭代過程，直到滿足每個變數的收斂條件，從而得到最終收斂的答案。

雅可比演算法

雅可比演算法如下：

• 從未知數（x）的一些初始猜測陣列開始。

• 透過在如下所示的方程重新排列形式中代入x的猜測值來計算新的x：

$$\mathrm{x_{i_{new}} \: = \: − \frac{1}{a_{i,i}}(\sum_{j=1,j \neq i}^n a_{i,j}x_{j_{guess}} \: − \: b_{i})}$$

現在，$\mathrm{x_{i_{new}}}$將是當前迭代中獲得的新x值。

• 下一步是評估新值和猜測值之間的誤差，即$\mathrm{\lvert x_{new} \: − \: x_{guess} \rvert}$。如果誤差大於某個收斂準則（我們將其設為$\mathrm{10^{−5}}$），則將新值賦給舊猜測值，即$\mathrm{x_{guess} \: = \: x_{new}}$，然後開始下一次迭代。

• 否則，$\mathrm{x_{new}}$是最終答案。

雅可比演算法——示例

讓我們透過以下示例來演示該演算法：

$$\mathrm{20x \: + \: y \: − \: 2z \: = \: 17}$$

$$\mathrm{3x \: + \: 20y \: − \: z \: = \: −18}$$

$$\mathrm{2x \: − \: 3y \: + \: 20z \: = \: 25}$$

將上述方程重新排列如下：

$$\mathrm{x_{new} \: = \: (−y_{guess} \: + \: 2z_{guess} \: + \: 17)/20}$$

$$\mathrm{y_{new} \: = \: (−3x_{guess} \: + \: z_{guess} \: − \: 18)/20}$$

$$\mathrm{z_{new} \: = \: (−2x_{guess} \: + \: 3y_{guess} \: + \: 25)/20}$$

現在，這些方程將在while迴圈中求解，以根據猜測值獲得未知數的新值。

實現雅可比方法的Python程式

下面顯示了實現雅可比方法（按方程實現）的程式：

示例

# Importing module for plotting and array
from pylab import *
from numpy import *

# Initial guess to start with
xg=0
yg=0
zg=0

# Setting error to move into the while loop
error=1

# Setting up iteration counter
count=0

while error>1.E-5:
   count+=1
    
   #Evaluating new values based on old guess
   x=(17-yg+2*zg)/20
   y=(zg-18-3*xg)/20
   z=(25-2*xg+3*yg)/20

   # Error evaluation and plotting
   error = abs(x-xg)+abs(y-yg)+abs(z-zg)
   figure(1,dpi=300)
   semilogy(count,error,'ko')
   xlabel('iterations')
   ylabel('error')

   # Updating the Guess for next iteration.
   xg=x
   yg=y
   zg=z

savefig('error_jacobi.jpg')
print(f'x={round(x,5)}, y={round(y,5)}, z={round(z,5)}')

輸出

程式輸出將是

$$\mathrm{x \: = \: 1.0 \: , \: y \: = \: −1.0 \: , \: z \: = \: 1.0}$$

下圖顯示了每次迭代次數的誤差變化：

可以觀察到，收斂解在第8次迭代時達到。

但是，我認為如果我們可以將線性方程組的元素輸入到矩陣中並以矩陣形式而不是方程形式求解它們，將會更好。因此，執行此任務的過程程式碼如下：

示例

# Importing module
from pylab import *
from numpy import *

#-----------------------------------------#
# Array of coefficients of x
a=array([[20,1,-2],[3,20,-1],[2,-3,20]])

# Array of RHS vector
b=array([17,-18,25])

# Number of rows and columns
n=len(b)
#-----------------------------------------#


# Setting up the initial guess array
xg=zeros(len(b))

# Starting error to enter in loop
error=1

# Setting iteration counter
count=0

# Generating array for new x
xn=empty(len(b)) 
#-----------------------------------------#


while error>1.E-5:
   count+=1
   for i in range(n):
      sum1=0
      for j in range(n):
         if i!=j:
            sum1=sum1+a[i,j]*xg[j]
      xn[i]=(-1/a[i,i])*(sum1-b[i])


   # Error evaluation and plotting
   error = sum(abs(xn-xg))
    
   figure(1,dpi=300)
   semilogy(count,error,'ko')
   xlabel('iterations')
   ylabel('error')

   # Substituting new value as the Guess for next iteration
   xg=xn.copy()
#-----------------------------------------#

savefig('error_jacobi.jpg')
print('x: ',xn)

輸出

x: [ 1.00000007 -0.99999983 1.00000005]

結論

在本教程中，我們解釋瞭如何使用Python來模擬雅可比迭代法求解聯立線性方程組。討論了兩種方法：方程法和矩陣法。