【數值計算方法】線性方程組迭代演算法的Python實現

FE-有限元鹰發表於2024-08-07
演算法Python

線性方程組迭代演算法的Python實現

jacobi,GS,SOR迭代法

def JacobiIter(A:np.ndarray,
                b:np.ndarray,
                tol:float=1e-5,
                maxIter:int=100)->Tuple[np.ndarray,np.ndarray]:
    """使用Jacobi迭代法求解線性方程組Ax=b
    
    input:
    A: np.ndarray, 係數矩陣
    b: np.ndarray, 右端常數
    tol: float, 誤差限
    maxIter: int, 最大迭代次數
    
    output:
    x: np.ndarray, 解向量
    errors: np.ndarray, 誤差序列
    """
    from numpy import dot
    from numpy.linalg import norm
    x0=np.zeros(b.shape)
    L=-1*np.tril(A,k=-1).copy()
    U=-1*np.triu(A,k=1).copy()
    D=np.diag(np.diag(A)).copy()
    Dinv=np.linalg.inv(D)
    errors=[]
    for i in range(maxIter):
        x_next=dot(Dinv,(dot((L+U),x0)+b))
        # error check
        error=norm(b-dot(A,x_next),2)/norm(b,2)
        errors.append(error)
        if error<tol:
            return x_next,np.array(errors)
        else:
            x0=x_next

def GaussIter(A:np.ndarray,
                b:np.ndarray,
                tol:float=1e-5,
                maxIter:int=100)->Tuple[np.ndarray,np.ndarray]:
    """使用Gauss-Seidel迭代法求解線性方程組Ax=b
    
    input:
    A: np.ndarray, 係數矩陣
    b: np.ndarray, 右端常數
    tol: float, 誤差限
    maxIter: int, 最大迭代次數
    
    output:
    x: np.ndarray, 解向量
    errors: np.ndarray, 誤差序列
    """
    from numpy import dot
    from numpy.linalg import norm
    x0=np.zeros(b.shape)
    L=-1*np.tril(A,k=-1).copy()
    U=-1*np.triu(A,k=1).copy()
    D=np.diag(np.diag(A)).copy()
    DsubLinv=np.linalg.inv(D-L)
    errors=[]
    for i in range(maxIter):
        x_next=DsubLinv.dot(U).dot(x0)+DsubLinv.dot(b)
        # error check
        error=norm(b-dot(A,x_next),2)/norm(b,2)
        errors.append(error)
        if error<tol:
            return x_next,np.array(errors)
        else:
            x0=x_next

def SORIter(A:np.ndarray,
                b:np.ndarray,
                w:float=1.0,
                tol:float=1e-5,
                maxIter:int=100)->Tuple[np.ndarray,np.ndarray]:
    """使用SOR迭代法求解線性方程組Ax=b
    
    input:
    A: np.ndarray, 係數矩陣
    b: np.ndarray, 右端常數
    w: float, 鬆弛因子(0~2.0)
    tol: float, 誤差限
    maxIter: int, 最大迭代次數
    
    output:
    x: np.ndarray, 解向量
    errors: np.ndarray, 誤差序列
    """
    from numpy import dot
    from numpy.linalg import norm
    x0=np.zeros(b.shape)
    L=-1*np.tril(A,k=-1).copy()
    U=-1*np.triu(A,k=1).copy()
    D=np.diag(np.diag(A)).copy()
    
    DsubOmegaLinv=np.linalg.inv(D-w*L)
    errors=[]
    for i in range(maxIter):
        x_next=DsubOmegaLinv.dot((1-w)*D+w*U).dot(x0)+w*DsubOmegaLinv.dot(b)
        # error check
        error=norm(b-dot(A,x_next),2)/norm(b,2)
        errors.append(error)
        if error<tol:
            return x_next,np.array(errors)
        else:
            x0=x_next
  • 驗證
import numpy as np
from formu_lib import *
A=np.array([[2,-1,0],
            [-1,3,-1],
            [0,-1,2]])
b=np.array([1,8,-5])
extractVal=np.array([2,3,-1])

x1,er1=JacobiIter(A,b)
x2,er2=GaussIter(A,b)
x3,er3=SORIter(A,b,1.2)

ind1,ind2,ind3=list(range(len(er1))),list(range(len(er2))),list(range(len(er3)))
plotLines([ind1,ind2,ind3],[er1,er2,er3],["Jacobi iter error","Gauss iter error","SOR iter error"])

showError(x1,extractVal)
showError(x2,extractVal)
showError(x3,extractVal)

img

# 雅可比迭代法
數值解: [ 1.9999746   2.99999435 -1.0000254 ],
精確解: [ 2  3 -1],
誤差: 9.719103983280175e-06
# GS迭代法
數值解: [ 1.9999619  2.9999746 -1.0000127],
精確解: [ 2  3 -1],
誤差: 1.2701315856479742e-05
# SOR迭代法
數值解: [ 2.00001461  2.999993   -1.00000098],
精確解: [ 2  3 -1],
誤差: 4.338862621105977e-06

正定對稱線性方程組的不定常迭代:最速下降法,共軛梯度法

def SPDmethodSolve(A:np.ndarray,
                    b:np.ndarray,
                    tol:float=1e-5,
                    maxIter:int=200)->Tuple[np.ndarray,np.ndarray]:
    """使用最速下降法求解線性方程組Ax=b
    
    input:
    A: np.ndarray, 係數矩陣,必須是對稱正定矩陣
    b: np.ndarray, 右端常數
    tol: float, 誤差限
    maxIter: int, 最大迭代次數
    
    output:
    x: np.ndarray, 解向量
    errors: np.ndarray, 誤差序列
    """
    from numpy import dot
    from numpy.linalg import norm
    x0=np.zeros(b.shape)
    i,errors=0,[]
    while True :
        if i>maxIter:
            maxIter=1.5*maxIter
            print(f"迭代次數過多,自動調整為 {maxIter}")
        # 計算殘量r^k作為前進方向.
        r=b-dot(A,x0)
        # 計算前進距離a_k
        a=InnerProduct(r,r)/InnerProduct(dot(A,r),r)
        x_next=x0+a*r
        error=norm(b-dot(A,x_next),2)/norm(b,2)
        errors.append(error)
        if errors[-1]<tol:
            return x_next,np.array(errors)
        else:
            x0=x_next
            i+=1
        
def conjGrad(A:np.ndarray,
                b:np.ndarray,
                tol:float=1e-5,
                maxIter:int=200)->Tuple[np.ndarray,np.ndarray]:
    """使用共軛梯度法求解線性方程組Ax=b
    
    input:
    A: np.ndarray, 係數矩陣,必須是對稱正定矩陣
    b: np.ndarray, 右端常數
    tol: float, 誤差限
    maxIter: int, 最大迭代次數
    
    output:
    x: np.ndarray, 解向量
    errors: np.ndarray, 誤差序列
    """
    from numpy import dot
    from numpy.linalg import norm
    # 選擇初值x0,初始方向p0=r0=b-Ax0
    x0=np.zeros(b.shape)
    i,errors=0,[]
    r0=b-dot(A,x0)
    p_0=b-dot(A,x0)
    errors.append(norm(r0,2)/norm(b,2))
    while True :
        if i>maxIter:
            maxIter=1.5*maxIter
            print(f"迭代次數過多,自動調整為 {maxIter}")
        # 計算a_k,x^{k+1}=x_k+a_k*p_k
        a_k=InnerProduct(r0,p_0)/InnerProduct(dot(A,p_0),p_0)
        x_next=x0+a_k*p_0
        # 計算下一步的殘量
        r_k_next=b-dot(A,x_next)
        errors.append(norm(r_k_next,ord=2)/norm(b,2))
        # 如果殘量足夠小,則停止迭代
        if errors[-1]<tol:
            return x_next,np.array(errors)
        else:
            # 計算下一步的搜尋方向
            beta_k=-1*InnerProduct(r_k_next,A.dot(p_0))/InnerProduct(p_0,A.dot(p_0))
            p_0=r_k_next+beta_k*p_0
            x0=x_next
            i+=1
  • 驗證

img

from formu_lib import *
import numpy as np

A=np.array([[4,-1,0],
            [-1,4,-1],
            [0,-1,4]])
b=np.array([3,2,3])
extractVal=np.array([1,1,1])

x1,er1=SPDmethodSolve(A,b,1e-6)
x2,er2=conjGrad(A,b,1e-6)

plotLines([list(range(len(er1))),list(range(len(er2)))],[er1,er2],["SPD method error","conjugate gradient error"])

showError(x1,extractVal)
showError(x2,extractVal)

img

# SPD method
數值解: [0.99999951 0.99999951 0.99999951],
精確解: [1 1 1],
誤差: 4.891480784863234e-07
# conjugate gradient method
數值解: [1. 1. 1.],
精確解: [1 1 1],
誤差: 0.0

相關文章