夜间模式暗黑模式
字体
阴影
滤镜
圆角
主题色
关于 Ax=b 的迭代法求解探究

实验内容

  1. 使用Python实现 (i)Jacobi (ii)Gauss-Seidel (iii)SOR 迭代法
  2. 随机生成(i)三角占优 (ii)非三角占优 矩阵
  3. 验证算法的收敛性
  4. 研究SOR中w对算法结果的影响

实验过程

Jacobi迭代法

根据如下推导过程:

\begin{aligned}
A x &=b \
(D+L+U) x &=b \
D x &=b-(L+U) x \
x &=D^{-1}(b-(L+U) x)
\end{aligned}

其中$D$为$A$中的主对角线部分,$L$为$A$中的下三角部分,$U$为$A$中的上三角部分

我们得到 Jacobi 迭代法的一般形式:

\begin{aligned}
x_{0} &=\text { initial vector } \
x_{k+1} &=D^{-1}\left(b-(L+U) x_{k}\right) \text { for } k=0,1,2, \ldots .
\end{aligned}

使用 python 代码实现如下:

#Jacobi method
def Jacobi(LplusU,inv_D,b,iteritor=100,detect=20):
    """'x0_list' records the value of x[0] in each iteration, val 'detect' checks the distance of adjacent iterations in 'detect' rounds """
    x = np.zeros((b.shape[0],1))
    x0_list = [0]
    for i in range (iteritor):
        x_old = x.copy()
        x = np.matmul(inv_D, b - np.matmul(LplusU,x))
        x0_list.append(x[0][0])
        if i%detect==0 and np.linalg.norm(x_old-x) < 0.0000001:
            break
    return x,x0_list

Gauss-Seidel 迭代法

Gauss-Seidel迭代法由Jacobi迭代法改进而来,其变化为对于每一轮迭代$i$,计算$x_i[k]$的过程中使用本轮迭代得到的 $ xi[j],j<k $ ,而非上一轮的结果 $ x{i-1}[j],j<k $ 。

其一般形式如下:

\begin{aligned}
x_{0} &=\text { initial vector } \
x_{k+1} &=D^{-1}\left(b-U x_{k}-L x_{k+1}\right) \text { for } k=0,1,2, \ldots
\end{aligned}

使用python代码实现如下:

#Gauss-Seidel method
def Gauss_Seidel(LplusU,inv_D,b,iteritor=100,detect=20):
    x = np.zeros((b.shape[0],1))
    x0_list = [0]
    for i in range(iteritor):
        x_old = x.copy()
        for j in range(x.shape[0]):
            x_tmp = np.matmul(inv_D, b - np.matmul(LplusU,x))
            x[j][0] = x_tmp[j][0]
        x0_list.append(x[0][0])
        if i%detect==0 and np.linalg.norm(x_old-x) < 0.0000001:
            break
    return x,x0_list

SOR(Successive Over-Relaxation)

SOR 由 Gauss-Seidel 迭代法变化而来,加入松弛参数$\omega$。设 Gauss-Seidel 迭代法在第$i$轮迭代中得到$x_i[k]$,则 SOR 迭代法在第$i$轮迭代中得到 $x’_i[k]=\omega xi[k]+(1-\omega)x{i-1}[k]$,即为本次与上次Gauss-Seidel迭代法的加权平均值。这样的变化引入了物理中“惯性” 的思想。

当松弛参数$\omega > 1$时被认为是超松弛(over-relaxation)的。

其一般形式如下:

\begin{aligned}
x_{0} &=\text { initial vector } \
x_{k+1} &=(\omega L+D)^{-1}\left[(1-\omega) D x_{k}-\omega U x_{k}\right]+\omega(D+\omega L)^{-1} b \text { for } k=0,1,2, \ldots .
\end{aligned}

使用python代码实现如下:

#SOR
def SOR(LplusU,inv_D,b,w=1,iteritor=100,detect=20):
    x = np.zeros((b.shape[0],1))
    x0_list = [0]
    for i in range(iteritor):
        x_old = x.copy()
        for j in range(x.shape[0]):
            x_tmp = np.matmul(inv_D, b - np.matmul(LplusU,x))
            x[j][0] = (1-w)*x[j][0]+w*x_tmp[j][0]
        x0_list.append(x[0][0])
        if i%detect==0 and np.linalg.norm(x_old-x) < 0.0000001:
            break
    return x,x0_list

生成随机矩阵

使用函数generate_matrix生成给定行数在给定范围内随机的,元素在给定范围内随机的,满足严格对角占优(非对角占优)的增广矩阵。

def isDiagDom(a):
    """Judge whether a matrix is diagonal dominant"""
    flag = True
    for i in range(a.shape[0]):
        absSum = 0
        for j in range(a.shape[0]):
            absSum += abs(a[i][j])
        if 2*abs(a[i][i]) <= absSum:
            flag = False
    return flag

def generate_matrix(row_lower_limit=2,row_upper_limit=20,element_lower_limit=-100,element_upper_limit=100,DiagDom=True):
    """generate a (not) diagonal dominant random matrix in a given range; 'count' records failed-generation counts"""
    if DiagDom:
        flag = True
        count = 0
        while flag:
            row = np.random.randint(row_lower_limit,row_upper_limit)
            a = np.random.randint(element_lower_limit,element_upper_limit,(row,row+1))
            if isDiagDom(a):
                flag = False
            else:
                count += 1
        return a,count
    else:
        flag = True
        count = 0
        while flag:
            row = np.random.randint(row_lower_limit,row_upper_limit)
            a = np.random.randint(element_lower_limit,element_upper_limit,(row,row+1))
            if not isDiagDom(a):
                flag = False
            else:
                count += 1
        return a,count

验证算法的收敛性

严格对角占优矩阵

我们已经知道以上迭代算法对于给定的$Ax=b$收敛的充分条件是矩阵$A$严格对角占优。因此我们首先使用随机生成的严格对角占优矩阵来测试各迭代算法的收敛性。

在此之前还要对随机生成的增广矩阵$a$进行预处理,得到迭代算法中需要的A的上三角与下三角部分LplusU,A的主对角线部分的逆inv_D​。

def preProcess(a):
    b = np.zeros((a.shape[0],1))
    for i in range(0,a.shape[0]):
        b[i][0] = a[i][a.shape[0]]
    A = np.resize(a.T,(a.shape[0],a.shape[0]))
    A = A.T
    D = np.zeros((A.shape[0],A.shape[0]))
    for i in range (0,D.shape[0]):
        D[i][i] = A[i][i]
    LplusU = A - D
    inv_D = np.linalg.inv(D)
    return A,b,LplusU,inv_D

接下来使用三种迭代算法验证算法的收敛性(这里SOR中$\omega=1.25$)

a,count = generate_matrix()
A,b,LplusU,inv_D = preProcess(a)
print('A:\n',A,'\n\n','b:\n',b,'\n\n','LplusU:\n',LplusU,'\n\n','inv_D\n',inv_D,'\n\n')

fig, ax = plt.subplots()

x,x0_list = Jacobi(LplusU,inv_D,b)
print("Jacobi iteration")
print("x:     iteration:",len(x0_list)-1,'\n',x,'\n')
ax.plot(range(len(x0_list)), x0_list,label='Jacobi')

print("Gauss_Seidel iteration")
x,x0_list = Gauss_Seidel(LplusU,inv_D,b)
print("x:     iteration:",len(x0_list)-1,'\n',x,'\n')
ax.plot(range(len(x0_list)), x0_list,label='Gauss_Seidel')

print("SOR iteration with w=1.25")
x,x0_list = SOR(LplusU,inv_D,b,w=1.25)
print("x:     iteration:",len(x0_list)-1,'\n',x,'\n')
ax.plot(range(len(x0_list)), x0_list,label='SOR with w=1.25')

print("Standard solution")
print(np.linalg.solve(A,b))

ax.legend()

得到一次结果如下:

A:
 [[ 74   4  33]
 [-56 -71  -4]
 [ 28 -37 -88]] 

 b:
 [[ -2.]
 [-97.]
 [ 73.]] 

 LplusU:
 [[  0.   4.  33.]
 [-56.   0.  -4.]
 [ 28. -37.   0.]] 

 inv_D
 [[ 0.01351351  0.          0.        ]
 [-0.         -0.01408451 -0.        ]
 [-0.         -0.         -0.01136364]] 


Jacobi iteration
x:     iteration: 41 
 [[ 0.42808818]
 [ 1.09351403]
 [-1.15310852]] 

Gauss_Seidel iteration
x:     iteration: 21 
 [[ 0.42808818]
 [ 1.09351403]
 [-1.15310852]] 

SOR iteration with w=1.25
x:     iteration: 100 
 [[ 0.4306864 ]
 [ 1.09040026]
 [-1.14948132]] 

Standard solution
[[ 0.42808818]
 [ 1.09351403]
 [-1.15310852]]

通过简单的分析不难得到,三种迭代算法均收敛于方程组的标准解处。且Gauss-Seidel,Jacobi,SOR的收敛速度依次下降。

经过多次测试,均得到三种算法收敛于解的结果,以此验证了算法在给定矩阵为严格对角占优时的收敛性。

非严格对角占优矩阵

下面使用随机生成的非严格对角占优矩阵来测试各迭代算法的收敛性。这里只需要对测试模块中generate_matrix的参数做出调整,故不再列出代码。

得到一次结果如下:

A:
 [[ -62  -74   82   73  -85   79  -85   73    8  -69   34]
 [  23   32   28  -95  -60   94   48  -33   58  -56   44]
 [  66   67  -91  -92  -41  -25  -50   66   40   70   19]
 [ -88  -64  -63   22   92  -25   38  -91 -100    8  -70]
 [ -68  -99  -68  -40  -46  -47  -99   55   16  -95   57]
 [ -29   78  -33   73  -56   62  -88   28   70  -81   95]
 [ -28  -12  -11  -69  -45   -3   66   63  -54   49   68]
 [ -29   18   82   21   71   66   98   -4    0    9  -54]
 [ -50   90  -97  -75   84  -37   32   19  -75   72   61]
 [  52    5   60   87   43  -89  -93  -85   60   44   32]
 [ -77   15  -84   25   37  -70  -99  -78  -22   10  -35]] 

 b:
 [[-78.]
 [-78.]
 [-36.]
 [ 63.]
 [ 21.]
 [ 94.]
 [ 14.]
 [ -8.]
 [ 62.]
 [ 48.]
 [-47.]] 

 LplusU:
 [[   0.  -74.   82.   73.  -85.   79.  -85.   73.    8.  -69.   34.]
 [  23.    0.   28.  -95.  -60.   94.   48.  -33.   58.  -56.   44.]
 [  66.   67.    0.  -92.  -41.  -25.  -50.   66.   40.   70.   19.]
 [ -88.  -64.  -63.    0.   92.  -25.   38.  -91. -100.    8.  -70.]
 [ -68.  -99.  -68.  -40.    0.  -47.  -99.   55.   16.  -95.   57.]
 [ -29.   78.  -33.   73.  -56.    0.  -88.   28.   70.  -81.   95.]
 [ -28.  -12.  -11.  -69.  -45.   -3.    0.   63.  -54.   49.   68.]
 [ -29.   18.   82.   21.   71.   66.   98.    0.    0.    9.  -54.]
 [ -50.   90.  -97.  -75.   84.  -37.   32.   19.    0.   72.   61.]
 [  52.    5.   60.   87.   43.  -89.  -93.  -85.   60.    0.   32.]
 [ -77.   15.  -84.   25.   37.  -70.  -99.  -78.  -22.   10.    0.]] 

 inv_D
 [[-0.01612903 -0.         -0.         -0.         -0.         -0.
  -0.         -0.         -0.         -0.         -0.        ]
 [ 0.          0.03125     0.          0.          0.          0.
   0.          0.          0.          0.          0.        ]
 [-0.         -0.         -0.01098901 -0.         -0.         -0.
  -0.         -0.         -0.         -0.         -0.        ]
 [ 0.          0.          0.          0.04545455  0.          0.
   0.          0.          0.          0.          0.        ]
 [-0.         -0.         -0.         -0.         -0.02173913 -0.
  -0.         -0.         -0.         -0.         -0.        ]
 [ 0.          0.          0.          0.          0.          0.01612903
   0.          0.          0.          0.          0.        ]
 [ 0.          0.          0.          0.          0.          0.
   0.01515152  0.          0.          0.          0.        ]
 [-0.         -0.         -0.         -0.         -0.         -0.
  -0.         -0.25       -0.         -0.         -0.        ]
 [-0.         -0.         -0.         -0.         -0.         -0.
  -0.         -0.         -0.01333333 -0.         -0.        ]
 [ 0.          0.          0.          0.          0.          0.
   0.          0.          0.          0.02272727  0.        ]
 [-0.         -0.         -0.         -0.         -0.         -0.
  -0.         -0.         -0.         -0.         -0.02857143]] 


Jacobi iteration
x:     iteration: 1000 
 [[nan]
 [nan]
 [nan]
 [nan]
 [nan]
 [nan]
 [nan]
 [nan]
 [nan]
 [nan]
 [nan]] 

Gauss_Seidel iteration
x:     iteration: 1000 
 [[nan]
 [nan]
 [nan]
 [nan]
 [nan]
 [nan]
 [nan]
 [nan]
 [nan]
 [nan]
 [nan]] 

SOR iteration with w=1.25
x:     iteration: 1000 
 [[nan]
 [nan]
 [nan]
 [nan]
 [nan]
 [nan]
 [nan]
 [nan]
 [nan]
 [nan]
 [nan]] 

Standard solution
[[ 0.3794794 ]
 [-0.75910988]
 [-1.31754506]
 [ 1.27187867]
 [ 0.72690897]
 [ 0.09068645]
 [ 1.03407261]
 [ 0.15584119]
 [ 0.74236809]
 [ 0.12988662]
 [ 1.13856391]]

通过简单的分析不难得到,三种迭代算法均发散。

经过多次测试,均得到三种算法发散的结果,以此验证了算法在给定矩阵为非严格对角占优时不具有收敛性。

总结

给定矩阵为严格对角占优矩阵,则三种算法收敛到方程组的解;反之三种算法均发散。

SOR中w对算法结果的影响

我们这里使用随机生成的严格对角占优矩阵对不同的$\omega$值进行测试,并分析结果。

a,count = generate_matrix()
A,b,LplusU,inv_D = preProcess(a)
print('A:\n',A,'\n\n','b:\n',b,'\n\n','LplusU:\n',LplusU,'\n\n','inv_D\n',inv_D,'\n\n')

fig, ax = plt.subplots()

print("SOR iteration")
for i in np.arange(0.5,1.5,0.1):
    x,x0_list = SOR(LplusU,inv_D,b,w=i,iteritor=100)
    print("x:     iteration:",len(x0_list)-1,'\n',x,'\n')
    ax.plot(range(len(x0_list)), x0_list,label='SOR with w='+str(i))

print("Standard solution")
print(np.linalg.solve(A,b))

ax.legend()

得到如下几组测试结果:




通过分析可以看出,当$\omega$取值为$[0.5,1.4]$中的某一个值$\omega_0$时SOR收敛速度最快,当$\omega<\omega_0$或$\omega>\omega_0$时收敛速度逐渐变慢,甚至出现发散的情况。此外,可能出现$\omega_0<1$的情况,即此时取$\omega>1$比$\omega=1$(Gauss-Seidel收敛)收敛速度慢,这表明使用SOR并不总是比使用Gauss-Seidel更加有效。

实验总结

该实验探究了(i)Jacobi (ii)Gauss-Seidel (iii)SOR 迭代法 在不同矩阵(是否为对角线占优)中的收敛性,并对SOR迭代法中$\omega$的不同取值对算法结果的影响进行了分析。

值得一提的是,实验代码generate_matrix函数中使用count变量记录成功随机生成一个对角占优矩阵前失败的次数。而count值从数十到数百不等,即随机生成一个矩阵,其为对角占优矩阵的概率在10%以下,而这个概率随着矩阵的阶的增大而急剧下降。在所有生成的对角占优矩阵中,有90%以上的均为2阶矩阵,剩下的均为3阶矩阵。未发现3阶以上的矩阵。

这个结果表明,对于一般的解方程组问题而言,上述迭代法可以起作用的概率不大,甚至可以说微乎其微。而对于特定的对角占优矩阵而言,迭代法有着巨大的优势。

暂无评论

发送评论 编辑评论


				
上一篇
下一篇