Python-线性最优解迭代方法

为什么80%的码农都做不了架构师？>>>   

定义矩阵乘法

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))

转载于:https://my.oschina.net/VenusV/blog/2874799