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 <stdio.h> #include <stdlib.h> #include <conio.h> #include <math.h>
#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]);
x_former[i]=x[i];
}
printf("\n");
for(i=1;i<=n;++i)
{
x[i]=A[i][n+1];
for(j=1;j<=n;++j)
if(j!=i) x[i]-=A[i][j]*x[j]; if(fabs(A[i][i])>PRECISION)
x[i]/=A[i][i];
else
return 1;
}
++k;
}while(IsSatisfyPricision(x,x_former,n) && k<MAX_Number);
if(k>=MAX_Number)
return 1;
else
{
printf("\nGauss-Seidel迭代次数为%d 次",k);
return 0;
}
}
int main() //主函数
{
int n;
float A[MAX_n][MAX_n],x[MAX_n];
printf("\n方阵维数n=");
scanf("%d",&n);
if(n>=MAX_n-1)
{
printf("\n\007n must <%d!",MAX_n);
exit(0);
}
MatrixInput(A,n,n+1);
if(Jacobi_(A,x,n))
printf("\nGauss-Seidel迭代失败!");
else
{
printf("\n结果:");
VectorOutput(x,n);
}
printf("\n\n\007Press any key to quit!\n");
getch();
}。