提示
程序只是对实际数学公式的一个反映,因此在下面我都先给数学公式再给函数,当然,函数是已经封装好的可以直接用的。

Jacobi迭代法 链接到标题

先来说说他的步骤

  1. 设定初始值x0
  2. 设定迭代次数k
  3. 设定精度eps
  4. 设定迭代次数max_iter
  5. 设定系数矩阵A和常数项列向量b

其数学迭代公式为 $$x^{(k+1)}i=\frac{\displaystyle b_i-\sum{j=1,j\neq i}^n a_{ij}x^{(k)}j}{a{ii}}$$

def jacobi(A, b, x0, max_iter=10):
    n = len(b)
    x = x0.copy()
    x_new = np.zeros_like(x0)
    history = [x0.copy()]
    for _ in range(max_iter):
        for i in range(n):
            s = np.dot(A[i, :], x) - A[i, i] * x[i]
            x_new[i] = (b[i] - s) / A[i, i]
        x = x_new.copy()
        history.append(x.copy())
    return x, np.array(history)
function x = jacobi(A, b, k)
    % A: 系数矩阵
    % b: 常数项列向量
    % k: 最大迭代次数
    n = length(b);
    x = zeros(n, 1);
    x_new = zeros(n, 1);
    
    for iter = 1:k
        for i = 1:n
            sum = b(i);
            for j = 1:n
                if i ~= j
                    sum = sum - A(i,j) * x(j);
                end
            end
            x_new(i) = sum / A(i,i);
        end
        x = x_new;
    end
end

Gauss-Seidel迭代法 链接到标题

其数学迭代公式为 $$x^{(k+1)}i=\frac{\displaystyle b_i-\sum{j=1,j\neq i}^n a_{ij}x^{(k+1)}j-\sum{j=1}^n a_{ij}x^{(k)}j}{a{ii}}$$

def gauss_seidel(A,b,x0,max_iter=10):
    n=len(b)#获取方程组的阶数
    x=x0.copy()#复制初始值
    x_new=np.zeros_like(x0)#初始化新的解
    history=[x0.copy()]#记录迭代过程
    for _ in range(max_iter):#这个带_的循环的意思是循环max_iter次,与不带_的循环不同,不需要用到循环变量
        for i in range(n):
            s=np.dot(A[i,:],x)-A[i,i]*x[i] #计算除了第i个元素之外的所有元素的乘积。
            x[i]=(b[i]-s)/A[i,i]#计算新的解
            #上面的两行代码其实是书本上端1公式2.5
        history.append(x.copy())#记录迭代过程
    return x,np.array(history)#返回最终解和迭代过程

SOR迭代法 链接到标题

SOR迭代法的数学迭代公式是Gauss-Seidel迭代法的一个变种,增加了一个松弛因子$\omega$,其数学迭代公式含$\omega$ 为 $$x^{(k+1)}i=(1-\omega)x^{(k)}i+\frac{\displaystyle \omega}{a{ii}}\left(b_i-\sum{j=1,j\neq i}^n a_{ij}x^{(k+1)}j-\sum{j=1}^n a_{ij}x^{(k)}_j\right)$$


def sor_solver(A, b, x_exact, omega, tol=5e-6, max_iterations=1000):
    n = len(b)
    x = np.zeros_like(b)
    iteration = 0
    iteration_process = []

    for iteration in range(max_iterations):
        x_new = np.copy(x)
        
        for i in range(n):
            sigma = 0
            for j in range(n):
                if j != i:
                    sigma += A[i, j] * x_new[j]
            x_new[i] = (1 - omega) * x[i] + (omega / A[i, i]) * (b[i] - sigma)
        
        iteration_process.append(np.copy(x_new))
        
        if np.linalg.norm(x_new - x_exact, ord=np.inf) < tol:
            break
        
        x = x_new

    return x, iteration + 1, iteration_process
function x = SOR(A, b, k, omega)
    % A: 系数矩阵
    % b: 常数项列向量
    % k: 最大迭代次数
    % omega: 松弛因子
    n = length(b);
    x = zeros(n, 1); 
    
    for iter = 1:k
        for i = 1:n
            sum1 = b(i);
            sum2 = 0;
            
            % 上三角部分
            for j = 1:i-1
                sum1 = sum1 - A(i,j) * x(j);
            end
            
            % 下三角部分
            for j = i+1:n
                sum2 = sum2 - A(i,j) * x(j);
            end
            

            x(i) = (1 - omega) * x(i) + (omega / A(i,i)) * (sum1 + sum2);
        end
    end
end