几种求解线性方程组的算法

一、高斯消元法

import numpy as np

def gauss(A, b):

    '''

    高斯消元法求解Ax=b，返回x向量

    '''

    m,n=A.shape[0],A.shape[1]

    B=np.concatenate((A,b.reshape(m,1)),axis=1)     # 构造增广矩阵

    for k in range(m):

        if B[k,k]==0:

            for i in range(k+1,m):

                if B[i,k]!=0:

                    B[[k,i],:]=B[[i,k],:]

        for j in range(k+1,n):

            B[j,:]=B[j,:]-(B[j,k]/B[k,k])*B[k,:]

    x = np.empty(n)

    for k in range(n-1,-1,-1):

        x[k]=(B[k,n]-B[k,k+1:n].dot(x[k+1:]))/B[k,k]

    return x

二、雅可比迭代法

关于几种迭代算法的迭代公式

 

def gauss_seidel(A, b):

    '''

    gauss_seidel迭代法求解Ax=b，返回x向量

    '''

    n=A.shape[0]

    x=np.zeros(n)

    res=np.linalg.norm(A.dot(x)-b)

    k=1

    while res>1e-5 and k<1000:

        for i in range(n):

            x[i]=x[i]+(b[i]-A[i,:].dot(x))/A[i,i]

            k=k+1

        res=np.linalg.norm(A.dot(x)-b)

    return x

三、SOR迭代法

def sor(A, b, w):

    '''

    SOR迭代法求解Ax=b，返回x向量

    '''

    n=A.shape[0]

    x=np.zeros(n)

    res=np.linalg.norm(A.dot(x)-b)

    k=1

    while res>1e-5 and k<1000:

        for i in range(n):

            x[i]=x[i]+w*(b[i]-A[i,:].dot(x))/A[i,i]

            k=k+1

        res=np.linalg.norm(A.dot(x)-b)

    return x

四、实例以及执行结果

A=np.array([[10,-1,-2],[-1,10,-2],[-1,-1,5]])

b=np.array([7.2,8.3,4.2])

print(gauss(A, b))

print(jacobi(A, b))

print(gauss_seidel(A, b))

print(sor(A, b, 0.5))