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

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)#返回最终解和迭代过程

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