Jacobi Gauss-Seidel迭代法

一、简介

考虑线性方程组Ax = b时，一般当A为低阶稠密矩阵时，用主元消去法解此方程组是有效方法。但是，对于由工程技术中产生的大型稀疏矩阵方程组（A的阶数很高，但零元素较多，例如求某些偏微分方程数值解所产生的线性方程组），利用迭代法求解此方程组就是合适的，在计算机内存和运算两方面，迭代法通常都可利用A中有大量零元素的特点。Jacobi、Gauss-Seidel迭代法就是解决其中问题的两种方法。
详解

二、实现

# -*- coding: utf-8 -*-
"""
Created on Fri Dec 16 10:49:27 2016

    Jacobi Gauss-Seidel

@author: Administrator
"""

from numpy import *

def jacobi(A,b,iter=50):
    m,n = A.shape
    if m != n:
        return
    x = zeros(n,dtype=float)
    temp_x =zeros(n,dtype=float)
    print 'Jacobi:'
    for i in range(0,iter):
        for j in range(0,n):
            temp_x[j] = x[j] + 1/A[j,j] * (b[j] - sum(A[j,:n] * x[:n]))
        print temp_x - x
        x[:n] = temp_x[:n]
    return x

def gauss_seidel(A,b,iter=50):
    m,n = A.shape
    if m != n:
        return
    x = zeros(n,dtype=float)
    temp_x =zeros(n,dtype=float)
    print 'Gauss-Seidel:'
    for i in range(0,iter):
        temp_x[:n] = x[:n] 
        for j in range(0,n):
            x[j] = x[j] + 1/A[j,j] * (b[j] - sum(A[j,:n] * x[:n]))
        print x -temp_x
    return x

if __name__ == '__main__':
    A = array([[5,-1,-3],[-1,2,4],[-3,4,15]],dtype=float)
    b = array([-2,1,10],dtype=float)
    jacobi(A,b,iter=100)
    gauss_seidel(A,b,iter=100)