Python实现Gauss-Seider迭代法（超全）

目录

1、代码实现

2、结果

---

1、代码实现

author:hewang
Email:207962168


import numpy as np
import matplotlib.pyplot as plt

def matvec(A,v):
    M,N = A.shape
    if N != M:
        print("Matrix is not square")
        return -1
    N1 = v.shape[0]
    if N1 != M:
        print("The dimension of matrix and vector is not reasonable ")
    result = np.zeros((N1,1),dtype = np.float32)
    for i in range(N):
        for j in range(N):
            result[i] += A[i,j]*v[j]
    return result
def product(v1,v2):
    result = 0
    M,N = v1.shape[0],v2.shape[0]
   
    if N != M:
        print("The dimension of vectors is not reasonable")
        exit(1)
    for i in range(M):
        result += v1[i]*v2[i]
    return result

def vecnorm2(v1):
    N = v1.shape[0]
    result = 0
    for i in range(N):
        result += v1[i]*v1[i]
    result = np.sqrt(result)
    return result
    

def cgsolver(A,b,x0,maxit:int,tol=1.0e-8):
    if tol < 0:
        print("argument tol is not good (negitive)")
    if maxit <= 0:
        maxit = 1000
    xk = x0
    rk = b-matvec(A,xk)
    dk = rk

    for i in range(maxit):
        print("iteration:{}, norm(x)= {}, norm(rk) = {}".format(i, vecnorm2(xk),vecnorm2(rk)))
        if vecnorm2(rk) <1.0e-16:
            print("A good approximation solution is obtained")
            return xk
        Adk = matvec(A,dk)
        alphak = product(rk,rk)/product(dk,Adk)

        xkp = xk + alphak*dk
        res = xkp-xk
        if vecnorm2(res) < tol:
            print("A good approximation solution is obtained")
            return xk 
        xk = xkp
        rkTrk = product(rk,rk)
        rk = rk -alphak*Adk
        betak = product(rk,rk)/rkTrk
        dk = rk + betak*dk
    print("Max iteration reachs")
    return xk

def cgsolver1(A,b,x0,maxit:int,tol=1.0e-8):
    if tol < 0:
        print("argument tol is not good (negitive)")
    if maxit <= 0:
        maxit = 1000
    xk = x0
    rk = b-np.matmul(A,xk)
    dk = rk

    for i in range(maxit):
        print("iteration:{}, norm(x)= {}, norm(rk) = {}".format(i, np.linalg.norm(xk),np.linalg.norm(rk)))
        if np.linalg.norm(rk) <1.0e-16:
            print("A good approximation solution is obtained")
            return xk
        Adk = np.matmul(A,dk)
        alphak = np.vdot(rk,rk)/np.vdot(dk,Adk)

        xkp = xk + alphak*dk
        res = xkp-xk
        if np.linalg.norm(res) < tol:
            print("A good approximation solution is obtained")
            return xk 
        xk = xkp
        rkTrk = np.vdot(rk,rk)
        rk = rk -alphak*Adk
        betak = np.vdot(rk,rk)/rkTrk
        dk = rk + betak*dk
    print("Max iteration reachs")
    return xk
def main():
    tol = 1.0e-8
    A = np.array([[2,2],[2,5]],dtype=np.float32)
    b = np.array([[6],[3]],dtype=np.float32)
    x0 = np.array([[1],[3]],dtype=np.float32)
    maxit = 20
    x= cgsolver1(A,b,x0,maxit,tol)
    print(x)




    

if __name__ == '__main__':
    main()

2、结果

iteration:0, norm(x)= 3.1622776985168457, norm(rk) = 14.142135620117188
iteration:1, norm(x)= 0.7820295691490173, norm(rk) = 3.8569462299346924
iteration:2, norm(x)= 4.123105525970459, norm(rk) = 2.454669584039948e-06
iteration:3, norm(x)= 4.123105525970459, norm(rk) = 8.890969382946423e-08
iteration:4, norm(x)= 4.123105525970459, norm(rk) = 6.228202090348411e-13
A good approximation solution is obtained
[[ 3.9999998]
 [-0.9999997]]