一、简介
考虑线性方程组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)