解线性方程组的迭代法_计算方法大作业

解线性方程组的迭代法

1.1 方法概述

对于线性方程组Ax=b ，其中A 为非奇异矩阵，当A 为大型稀疏矩阵时，考虑用迭代法求解上述方程组，其基本思想是求不动点

f BX X k k +=+)()1(, 即构造一个向量系列(){}

n X ，使其收敛至某个极限*X ，

则*X 就是要求的方程组的准确解。这儿主要介绍Jacobi 迭代法和Gauss-Seidel 迭代法的算法理论及数值实验。 1）Jacobi 迭代法算法概述

Jacobi 迭代法推导过程为将方程组：

⎪⎪⎩⎪⎪⎨⎧=+++=+++=+++n

n nn n n n n n n b x a x a x a b x a x a x a b x a x a x a 22112

22221211

1212111 )1(

在假设0≠ii a 下，改写成

()⎪⎪

⎩

⎪⎪⎨

⎧++++=++++=++++=--n n n n n n n n n n n g x b x b x b x g x b x b x b x g x b x b x b x 1122112

23231212113132121 )2( 如果引用系数矩阵

⎥⎥

⎥⎦

⎤

⎢⎢⎢⎣⎡=nn n n a a a a A 1111, ⎥⎥

⎥⎦

⎤

⎢⎢⎢⎣⎡=0011 n n b b B

及向量 ⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=n x x X 1，⎥⎥⎥

⎦

⎤⎢⎢⎢⎣⎡=n b b b 1，⎥⎥⎥⎦⎤

⎢⎢⎢⎣⎡=n g g g 1，

方程组（1）和（2）分别可写为：b AX =及g BX X +=，这样得到了jacobi

迭代格式01

f BX X

k k +=+用jacobi 迭代解方程组b AX =时，就可任意取

初值0X 带入迭代可知式g BX X k k +=+1，然后求k

k X ∞

→lim 。但是，n 比较大的时候，写方程组)1(和)2(是很麻烦的，如果直接由A ，b 能直接得到B ，g 就是矩阵与向量的运算了，那么如何得到B ，g 呢

实际上，如果引进非奇异对角矩阵()0≠ii a

⎥⎥

⎥⎦⎤

⎢⎢⎢⎣⎡=nn a a D 00011

将A 分解成：,D D A A +-=要求b AX =的解，实质上就有

,)(DX X D A AX +-=而D 奇异的，所以1-D 存在，

,)(X A D AX DX -+=从而有,11b D AX D X --+=我们在这里不妨令,1A D I B --=b D g 1-=就得到jacobi 迭代格式：f BX X k k +=+1 。 用向量的分量来表示为：

⎩

⎨⎧[]

,...

,,k ,

n ,...,i x a b

a x

n

i

j j )

k (j j i i

ii

)k (i

21021111==∑-=≠=+

其中()

()

()

()

()T

n

x ,...x ,x x

002010=为初始向量.

由此看出,雅可比迭代法公式简单,每迭代一次只需计算一次矩阵和向量的乘法.在电算时需要两组存储单元,以存放()

k x 及()

1+k x .

2）Gauss-Seidel 迭代法算法概述

由雅可比迭代公式可知,在迭代的每一步计算过程中是用()k x 的全部分量来计算()1+k x 的所有分量,显然在计算第i 个分量()

1+k i x 时,已经

计算出的最新分量()

()

11

11

+-+k i k x ,...,x 没有被利用，从直观上看,最新计算

出的分量可能比旧的分量要好些.因此，对这些最新计算出来的第

1+k 次近似()

1+k x

的分量()

1+k j

x 加以利用,就得到所谓解方程组的高斯—

塞德（Gauss-Seidel ）迭代法。

把矩阵A 分解成 U L D A --= 其中()nn

a ,...,a ,a diag D 2211=,U ,L --分别为A 的主对角元除外的下

三角和上三角部分,于是,方程组(1)便可以写成 ()b Ux x L D +=-

即 22f x B x += 其中 ()()b L D f ,

U L D B 1

21

2---=-=

以2B 为迭代矩阵构成的迭代法(公式)

()()221f x B x k k +=+

称为高斯—塞德尔迭代法(公式),用 量表示的形式为

⎩

⎨⎧[],...

,,k ,n ,,i x a x a b a x i j n i j )k (j ij )

k (j ij i ii )k (i

21021111111==∑∑--=-=+=++ 由此看出,高斯—塞德尔迭代法的一个明显的优点是,在电算时,只需一组存储单元(计算出()

1+k i x 后()k i x 不再使用,所以用()

1+k i

x 冲掉()

k i x ,以

便存放近似解. 1.2算法程序

1）Jacobi 迭代法

现在考虑Jacobi 迭代法的计算程序，按照算法(Jacobi 迭代法)编写Matlab 程序(Jacobi.m)如下：

function [x, k, index]=Jacobi(A, b, ep, it_max)

%求解线性方程组的Jacobi 迭代法，其中 % A ---方程组的系数矩阵 % b ---方程组的右端项

% ep ---精度要求。省缺为1e-5

% it_max ---最大迭代次数，省缺为100 % x ---方程组的解 % k ---迭代次数

% index --- index=1表示迭代收敛到指定要求；