Gauss-Seidel迭代法求解线性方程组

Gauss-Seidel迭代法求解线性方程组

一． 问题描述

用Gauss-Seidel 迭代法求解线性方程组

由Jacobi 迭代法中，每一次的迭代只用到前一次的迭代值。使用了两倍的存储空间，浪费了存储空间。若每一次迭代充分利用当前最新的迭代值，即在计算第i 个分量

)

1(+k i

x 时，用最新分量

)

1(1

+k x ，

⋅⋅⋅+)

1(2

k x )

1(1

-+k i x 代替旧分量

)

(1

k x ，

⋅⋅⋅)

(2

k x )

(1

-k i x ，可以起

到节省存储空间的作用。这样就得到所谓解方程组的Gauss-Seidel 迭代法。

二． 算法设计

将A 分解成U D L A --=，则b x =A 等价于b x =--U)D (L 则Gauss-Seidel 迭代过程

)

()1()1(k k k Ux Lx b Dx ++=++

故

)

()1()(k k Ux b x L D +=-+ 若设1

)(--L D 存在，则

b

L D Ux L D x k k 1)(1)1()()(--+-+-= 令

b

L D f U L D G 11)()(---=-=，

则Gauss-Seidel 迭代公式的矩阵形式为

f Gx x k k +=+)

()

1(

其迭代格式为

T

n x x x x )

()0()0(2)0(1)0(，，，⋅⋅⋅= (初始向量)，

)

(1

1

1

1

1

)()1()1(∑∑-=-+=++--=i j i i j k j ij

k j ij i ii i i

x a

x a b a x

)210i 210(n k ⋅⋅⋅=⋅⋅⋅=，，，；，，，

或者

⎪⎩

⎪⎨⎧--=⋅⋅⋅=⋅⋅⋅==∆+=∑∑-=-+=+++)

(1)210i 210(111

1)()1()1()()1(i j i i j k j ij k j ij i ii i i i k i k i x a x a b a x n k k x x x ，，，；，，，

三． 程序框图

四． 结果显示

TestBench

利用Gauss-Seidel 迭代法求解下列方程组

⎪⎩⎪

⎨⎧=+-=++--=++3

103220241225321

321321x x x x x x x x x ， 其中取→

=1)

0(x

。

运行程序

开始

读入数据n ，初始向量，增广矩阵

k

⇒1

i

ii i j i i j k j ij k j

ij i y a x a x

a b ⇒--∑∑-=-+=+/)(1

1

1

1

)()

1(

?

max ε<-i i y x

≥

k=N?

n

i x y k k i

i ，，⋅⋅⋅=⇒⇒+,3211

=

输出迭代失败标志

结束

输出

n

y y y ⋅⋅⋅21,

<

≠

依次输入：

1．方阵维数

2．增广矩阵系数

3．初始向量

得到：

迭代12次后算出x[1] = -4.0

x[2] = 3.0

x[3] = 2.0

五． 程序

#include #include #include #include

#define MAX_n

100

#define PRECISION 0.0000001 #define MAX_Number 1000

void VectorInput(float x[],int n) //输入初始向量 { int i;

for (i=1;i<=n;++i) { printf("x[%d]=",i); scanf("%f",&x[i]);

}

}

void MatrixInput(float A[][MAX_n],int m,int n) //输入增广矩阵

{

int i, j;

printf("\n输入系数矩阵：\n");

for(i=1;i<=m;++i)

{

printf("增广矩阵行数%d : ",i);

for(j=1;j<=n;++j)

scanf("%f",&A[i][j]);

}

}

void VectorOutput(float x[],int n) //输出向量

{

int i;

for(i=1;i<=n;++i)

printf("\nx[%d]=%f",i,x[i]);

}

int IsSatisfyPricision(float x1[],float x2[],int n) //判断是否在规定精度内{

int i;

for(i=1;i<=n;++i)

if(fabs(x1[i]-x2[i])>PRECISION) return 1;

return 0;

}

int Jacobi_(float A[][MAX_n],float x[],int n) //具体计算

{

float x_former[MAX_n];

int i,j,k;

printf("\n初始向量x0:\n");

VectorInput(x,n);

k=0;

do{

for(i=1;i<=n;++i)

{

printf("\nx[%d]=%f",i,x[i]);