algorithm-in-python/math/numericalAnalysis/solve-linear-by-iteration.py

''' mbinary
#########################################################################
# File : solve-linear-by-iteration.py
# Author: mbinary
# Mail: zhuheqin1@gmail.com
# Blog: https://mbinary.xyz
# Github: https://github.com/mbinary
# Created Time: 2018-10-02  21:14
# Description:
#########################################################################
'''

'''
#########################################################################
# File : solve-linear-by-iteration.py
# Author: mbinary
# Mail: zhuheqin1@gmail.com
# Blog: https://mbinary.xyz
# Github: https://github.com/mbinary
# Created Time: 2018-05-04  07:42
# Description: 
#########################################################################
'''
import numpy as np 
from operator import le,lt

def jacob(A,b,x,accuracy=None,times=6):
    ''' Ax=b,  arg x is the init val, times is the time of iterating'''
    A,b,x = np.matrix(A),np.matrix(b),np.matrix(x)
    n,m = A.shape 
    if n!=m:raise Exception("Not square matrix: {A}".format(A=A))
    if b.shape !=( n,1) : raise Exception('Error: {b} must be {n} x1 in dimension'.format(b = b,n=n))
    D = np.diag(np.diag(A))
    DI = np.zeros([n,n])
    for i in range(n):DI[i,i]= 1/D[i,i]
    R = np.eye(n) - DI * A
    g = DI * b
    print('R =\n{}'.format(R))
    print('g =\n{}'.format(g))
    last = -x
    if accuracy != None:
        ct=0
        while 1:
            ct+=1
            tmp = x-last
            last = x
            mx = max ( abs(i) for i in tmp) 
            if mx<accuracy:return x
            x = R*x+g
            print('x{ct} =\n{x}'.format(ct = ct,x=x))
    else:
        for i in range(times):
            x = R*x+g
            print('x{ct} =  \n{x}'.format(ct=i+1,x=x))
    print('isLimitd: {}'.format(isLimited(A)))
    return x
def gauss_seidel(A,b,x,accuracy=None,times=6):
    ''' Ax=b,  arg x is the init val, times is the time of iterating'''
    A,b,x = np.matrix(A),np.matrix(b),np.matrix(x)
    n,m = A.shape 
    if n!=m:raise Exception("Not square matrix: {A}".format(A=A))
    if b.shape !=( n,1) : raise Exception('Error: {b} must be {n} x1 in dimension'.format(b = b,n=n))
    D =np. matrix(np.diag(np.diag(A)))
    L = np.tril(A) - D  # L = np.triu(D.T) - D
    U = np.triu(A) - D
    DLI = (D+L).I
    S = - (DLI) * U
    f = (DLI)*b
    print('S =\n{}'.format(S))
    print('f =\n{}'.format(f))
    last = -x
    if accuracy != None:
        ct=0
        while 1:
            ct+=1
            tmp = x-last
            last = x
            mx = max ( abs(i) for i in tmp) 
            if mx<accuracy:return x
            x = S*x+f
            print('x{ct} =\n{x}'.format(ct=ct,x=x))
    else:
        for i in range(times):
            x = S*x+f
            print('x{ct} =  \n{x}'.format(ct=i+1,x=x))
    print('isLimitd: {}'.format(isLimited(A)))
    return x


def isLimited(A,strict=False):
    '''通过检查A是否是[严格]对角优来判断迭代是否收敛, 即对角线上的值是否都大于对应行(或者列)的值'''
    diag = np.diag(A)
    op = lt if strict else le
    if op(A.max(axis=0),diag).all(): return True 
    if op(A.max(axis=1), diag).all(): return True 
    return False

testcase=[]
def test():
    for func,A,b,x,*args in testcase:
        acc =None 
        times = 6
        if args !=[] :
            if isinstance(args[0],int):times = args[0]
            else : acc = args[0]
        return func(A,b,x,acc,times)


if __name__ =='__main__':
    A = [[2,-1,-1],
         [1,5,-1],
         [1,1,10]
        ]
    b = [[-5],[8],[11]]
    x = [[1],[1],[1]]
    #testcase.append([gauss_seidel,A,b,x])

    A = [[2,-1,1],[3,3,9],[3,3,5]]
    b = [[-1],[0],[4]]
    x = [[0],[0],[0]]
    #testcase.append([jacob,A,b,x])

    A = [[5,-1,-1],
         [3,6,2],
         [1,-1,2]
        ]
    b=  [[16],[11],[-2]]
    x = [[1],[1],[-1]]
    testcase.append([gauss_seidel,A,b,x,0.001])
    test()