Python-線性最優解迭代方法

定義矩陣乘法

def mult(h, x):
    result = []
    for x in h:
        summ = 0
        for index, y in enumerate(x):
            summ += y * x[index]
        result.append(summ)
    return result

建立希爾伯特矩陣

def create_hobert(n):
    h=[]
    for x in range(1, n + 1):
        row = []
        for y in range(1, n + 1):
            row.append(1 / (x + y - 1))
        h.append(row)
    return h

雅克比迭代法

def jacobi(A,B,sigma,N):
    n = len(A)
    x0 = []
    x = []
    for i in range(0,n):
        x0.append(0)
        x.append(0)
    for k in range(1,N+1):
        R = 0
        for i in range(0,n):
            sum_ax = 0
            for j in range(0,n):
                sum_ax = sum_ax + A[i][j] * x0[j]
            x[i] = x0[i] + ((B[i] - sum_ax)/A[i][i])
            if abs(x[i] - x0[i]) > R:
                R = abs(x[i] - x0[i])
        if R <= sigma:
            print("精確度等於",sigma,"時，jacobi迭代法須要迭代",k,"次收斂")
            return (x,k)
        for i in range(0,n):
            x0[i] = x[i]
    return (x,k)

Hauss-Seidel迭代法

def gauss_seidel(A,B,sigma,N):
    n = len(A)
    x0 = []
    x = []
    for i in range(0,n):
        x0.append(0)
        x.append(0)
    for k in range(1,N+1):
        R = 0
        for i in range(0,n):
            sum_ax = 0
            for j in range(0,n):
                if j >= i:
                    sum_ax = sum_ax + A[i][j] * x0[j]
                else:
                    sum_ax = sum_ax + A[i][j] * x[j]
            x[i] = x0[i] + ((B[i] - sum_ax)/A[i][i])
            if abs(x[i] - x0[i]) > R:
                R = abs(x[i] - x0[i])
        if R <= sigma:
            print("精確度等於",sigma,"時，gauss_seidel迭代法須要迭代",k,"次收斂")
            return (x,k)
        for i in range(0,n):
            x0[i] = x[i]
    return (x,k)

SOR迭代法

def sor(A,B,sigma,N,omega):
    n = len(A)
    x0 = []
    x = []
    for i in range(0,n):
        x0.append(0)
        x.append(0)
    for k in range(1,N+1):
        R = 0
        for i in range(0,n):
            sum_ax = 0
            for j in range(0,n):
                if j >= i:
                    sum_ax = sum_ax + A[i][j] * x0[j]
                else:
                    sum_ax = sum_ax + A[i][j] * x[j]
            x[i] = x0[i] + omega * ((B[i] - sum_ax)/A[i][i])
            if abs(x[i] - x0[i]) > R:
                R = abs(x[i] - x0[i])
        if R <= sigma:
            print("精確度等於",sigma,"時，鬆弛因子爲",omega,"時,sor迭代法須要迭代",k,"次收斂")
            print(x)
            return (x,k)
        for i in range(0,n):
            x0[i] = x[i]
    return (x,k)

測試希爾伯特矩陣H

if __name__ == "__main__":
    n = 10
    H = create_hobert(n)
    x = [1 for x in range(n)]
    b = mult(H,x)
    print('雅克比迭代法:',jacobi(H, b, 0.1, 200))
    print('SOR迭代法:',sor(H, b, 0.001, 20000, 1.5))
    # print('Guass-Seidel迭代法:',gauss_seidel(H,b,0.00001,20))