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

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

jacobi迭代法

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

• 驗證

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

x,er=JacobiIter(A,b)
plotLines([list(range(len(er))),],[er,],["Jacobi iter error"])

showError(x,extractVal)

數值解: [ 1.9999746   2.99999435 -1.0000254 ],
精確解: [ 2  3 -1],
誤差: 9.719103983280175e-06