美文网首页
只要一杯秋天的奶茶,就能学会Python数值分析(2)

只要一杯秋天的奶茶,就能学会Python数值分析(2)

作者: zengw20 | 来源:发表于2020-09-25 22:21 被阅读0次

    上节(https://www.jianshu.com/p/671a94ce586b)说到高斯消元法,今天从高斯列主元消元法开始,拓展到线性方程组的两种迭代方法:雅可比迭代法和高斯-赛德尔迭代法。同样,能力有限,希望读者指出我的问题,用代码和公式和我深入交流。

    2.列主元消去法

    参考的教材是《数值分析》(李庆扬等)的第148页到151页。
    列主元消去法是使用于主元为0或很小的情况,思路是选取这一列中绝对值最大的值所在的行与当前行进行替换,具体可以参考教材。再写代码方面,有了第一节的东西,写列主元消去法就很简单了,只要在高斯消去法之前加上换行就行。代码如下:

    def pivotTriangularMatrix(matA):
        '''
        先找主元,找到后跟高斯顺序消元操作一样
        '''
        numRows = matA.shape[0]
        for row in np.arange(0,numRows-1):
            Index = np.argmax(abs(matA[:,row]))
            matA[[row,Index],:] = matA[[Index,row],:] # 交换。也可以用np.copy来做交换
            
            # 交换行结束,开始消元
            for k in range(row+1,numRows):  # 在第row行的情况下,遍历剩下的行
                if matA[k,row] != 0:
                #第k(k小于row)行的新值等于该行减去第row行的值乘于一个系数。这个系数存在目的就是消元
                    matA[k,:]=matA[k,:]-(matA[k,row]/matA[row,row])*matA[row,:]
        return matA # 此时输入的矩阵也被改变  
    

    用教材148页的例题4来测试一下这个代码。

    test_mat3 = np.array([[0.001,2.0,3.0,1.0],[-1.0,3.712,4.623,2.0],[-2.0,1.072,5.643,3.0]],dtype=float)
    result = Gauss_solve(pivotTriangularMatrix(test_mat3)) # Gauss_solve函数是上一节的函数
    # 得出结果
    result = array([[-0.49039646],
                    [-0.05103518],
                    [ 0.36752025]])
    

    得出的结果和书中给出的答案是吻合的。暂时没有找其它算例测试。

    3.今日加餐

    今天舍友问我下面这个代码。刚好今天没去旁听。主要问这里面提到的独特的索引方式是什么。其实就是一组布尔值。这里举个例子,具体解释在代码里面。


    实现取斜对角线.png
    i,n = 1,3
    idx1 = np.arange(i,n+i)
    '''
    假设i=1,n=3
    idx1 = array([1, 2, 3])
    '''
    idx1>n-1
    '''
    输出为 array([False, False,  True]),这样就拿到了True所在的索引2
    那么idx1[idx1>n-1]在这里就相当于idx1[2],这个值等于3
    '''
    # 举例
    test_mat1 = np.array([[2,0,6],[1,4,3],[3,2,1]],dtype=float)
    test_mat[idx0,idx1]
    '''
    测试矩阵为:
    test_mat1 = array([[2., 0., 6.],
                       [1., 4., 3.],
                       [3., 2., 1.]])
    输出值为 array([0., 3., 3.]) 。实现了上图中取斜对角的操作
    '''
    

    二、线性方程组的迭代求解法

    1.雅可比迭代法

    该处参考的是《数值分析》(李庆扬等)第6章的6.1.1和6.2.1。书上关于雅可比迭代的内容有更多细节。这里我按照我对Gradient Descent的理解来编写该处代码。
    这是一个迭代的计算方法,从我的理解来看,这类方法与线性方程组直接求解法有2个较大区别:1.对于解,会有一个初始值,会从这个初始值开始按照迭代公式运算下去。2.迭代法都要有迭代结束条件,要么是迭代解与真实解之间的差小于人为设定的一个阈值,要么是达到人为设定的最大迭代次数。
    以下是我的代码,参考的是《数值分析》(李庆扬等)188页2.3这个公式,这个公式是一个展开的公式,我代码里写成矩阵运算的形式,这样可以减少使用for循环。

    
    def Jocobi(A,b,initial,delta,maxTimes):
        '''
        输入:A是系数矩阵,N阶方阵
              b是N*1列向量
              initial是解的初始值,N*1大小
              delta是人为设置的一个阈值,指到后期两次迭代之间解之间的差距已经很小了,用作迭代结束条件
              maxTimes指人为设定的最大迭代次数
              
              教材6.2.1将矩阵A分解成D,L,U,这里与之对应
        输出:迭代后的解析解
        '''
        D = np.diag(np.diag(A)) # 获得D矩阵
        L = np.tril(A,-1)       # 获得L矩阵
        U = np.triu(A,1)        # 获得U矩阵
        d = np.linalg.inv(D)    # 对D矩阵求逆
        G = -np.dot(d,L+U)      # D的逆矩阵乘于(L+U)矩阵
        f = np.dot(d,b)         # D的逆矩阵b向量
        X = np.dot(G,initial)+f # 初次的解
        for i in range(maxTimes):
            if np.linalg.norm(X - initial) > delta:
                initial = X
                X = np.dot(G,initial) + f
        return X
    

    利用第180页的例题1来对这个函数进行验证。解的初始值设为0解。

    import numpy as np
    A = np.array([[8,-3,2],[4,11,-1],[6,3,12]],dtype=float)
    b = np.array([[20],[33],[36]],dtype=float)
    initial = np.zeros((3,1))
    '''
    A = array([[ 8., -3.,  2.],
               [ 4., 11., -1.],
               [ 6.,  3., 12.]])
    b = array([[20.],
               [33.],
               [36.]])
    initial = array([[0.],
                    [0.],
                    [0.]])
    '''
    # 迭代求解
    X = Jocobi(A,b,initial,0.001,10)
    '''
    X = array([[3.00028157],
               [1.99991182],
               [0.99974048]])
    '''
    

    进过10次迭代后解出的结果与书上181页的结果一致。

    2.高斯-赛德尔迭代法

    有了雅可比迭代的代码,写高斯赛德尔迭代代码就非常简单,你可能会惊讶,这两个代码怎么是一样的。看下图可以看到这两的区别。参考教材的6.2.2节。


    image.png

    以下是我的代码:

    def Seidel(A,b,initial,delta,maxTimes):
        '''
        输入:A是系数矩阵,N阶方阵
              b是N*1列向量
              initial是解的初始值,N*1大小
              delta是人为设置的一个阈值,指到后期两次迭代之间解之间的差距已经很小了,用作迭代结束条件
              maxTimes指人为设定的最大迭代次数
              
              教材6.2.1将矩阵A分解成D,L,U,这里与之对应
        输出:迭代后的解析解
        '''
        # 先分解
        D = np.diag(np.diag(A)) 
        L = np.tril(A,-1)
        U = np.triu(A,1)
        
        d = np.linalg.inv( D + L ) # 这里是和雅克比不同的地方
        
        G = -np.dot(d,U)
        f = np.dot(d,b)
        X = np.dot( G ,initial ) + f
        for i in range(maxTimes):
            if np.linalg.norm( X - initial ) > delta:
                initial = X
                X = np.dot( G,initial )+f       
        return X
    

    同样用前面雅可比的算例来测试一下。

    X = Seidel(A,b,initial,0.0001,10)
    '''
    X = array([[3.00000201],
               [1.9999987 ],
               [0.99999932]])
    '''
    

    同样的例题,在同样迭代10次的情况下,与前面雅可比迭代的结果有一点点偏差,但这都在正常范围之内。与书上的结果基本吻合。

    今日小结

    今天可太忙了,从早到晚。时间全靠课间挤,关于迭代法也没有考虑收敛条件的问题。这个系列对我来说还挺好玩,打算有始有终做下去。原来,写东西(还不用负责)这件事情能这么快乐。
    今晚毕,明天继续。。。。。

    相关文章

      网友评论

          本文标题:只要一杯秋天的奶茶,就能学会Python数值分析(2)

          本文链接:https://www.haomeiwen.com/subject/deaxuktx.html