一. 问题描述
用Gauss-Seidel迭代法求解线性方程组
由Jacobi迭代法中,每一次的迭代只用到前一次的迭代值。使用了两倍的存储空间,浪
x费了存储空间。若每一次迭代充分利用当前最新的迭代值,即在计算第i个分量i用最新分量x1(k?1)(k?1)时,
,x2(k?1)???xi-1(k?1)代替旧分量x1,x2(k)(k)???xi-1,可以起到节省存储
(k)空间的作用。这样就得到所谓解方程组的Gauss-Seidel迭代法。
二. 算法设计
将A分解成A?L?D?U,则Ax?b等价于(L?D?U)x?b 则Gauss-Seidel迭代过程
Dx(k?1)?b?Lx(k?1)?Ux(k)故
(D?L)x(k?1)?b?Ux(k)
若设(D?L)存在,则
?1x(k?1)?(D?L)?1Ux(k)?(D?L)?1b
令
G?(D?L)?1U,f?(D?L)?1b
则Gauss-Seidel迭代公式的矩阵形式为
x(k?1)?Gx(k)?f
其迭代格式为
(0)(0)Tx(0)?(x1(0),x2,???,xn) (初始向量),
x或者
(i?1)ii?1i?11(k?1)1,2,???;i?0,1,2,???n) ?(bi??aijxj??aijx(jk)) (k?0,aiij?1j?i?1?xi(k?1)?xi(k)??xi(k?k?0,1,2,???;i?0,1,2,???n)?i?1i?11 ? (i?1)(k?1)xi?(bi??aijxj??aijx(jk))?aiij?1j?i?1?
三. 程序框图
开始 读入数据n,初始向量,增广矩阵 1?k(bi??axj?1i?1(k?1)ijj??aijx(jk))/aii?yij?i?1i?1 maxxi?yi???k?1?kyi?xii?1,2,3???,n?? k=N? ??输出输出迭代失败标志 y1,y2???yn 结束
四. 结果显示
TestBench
利用Gauss-Seidel迭代法求解下列方程组
?5x1?2x2?x3??12??(0)??x1?4x2?2x3?20, 其中取x?1。 ?2x?3x?10x?323?1运行程序
依次输入:
1.方阵维数
2.增广矩阵系数 3.初始向量
得到:
迭代12次后算出 x[1] = -4.0 x[2] = 3.0 x[3] = 2.0
五. 程序
#include
#define MAX_n
100
#define PRECISION 0.0000001 #define MAX_Number 1000
void VectorInput(float x[],int n) //输入初始向量 { }
void MatrixInput(float A[][MAX_n],int m,int n) //输入增广矩阵 { }
void VectorOutput(float x[],int n) //输出向量 { }
int IsSatisfyPricision(float x1[],float x2[],int n) //判断是否在规定精度内 {
int i; int i;
for(i=1;i<=n;++i)
printf(\,i,x[i]); int i, j;
printf(\输入系数矩阵:\\n\); for(i=1;i<=m;++i) { }
printf(\增广矩阵行数%d : \,i); for(j=1;j<=n;++j)
scanf(\,&A[i][j]);
for(i=1;i<=n;++i) { }
printf(\,i); scanf(\,&x[i]); int i;
}
int Jacobi_(float A[][MAX_n],float x[],int n) //具体计算 { do{ }
int main() //主函数 {
if(k>=MAX_Number) { }
printf(\迭代次数为%d 次\,k); return 0; return 1; else
for(i=1;i<=n;++i) { }
printf(\); for(i=1;i<=n;++i) { } ++k;
x[i]=A[i][n+1]; for(j=1;j<=n;++j)
if(j!=i) x[i]-=A[i][j]*x[j]; x[i]/=A[i][i]; return 1;
if(fabs(A[i][i])>PRECISION)
printf(\,i,x[i]); x_former[i]=x[i];
k=0;
printf(\初始向量x0:\\n\); VectorInput(x,n); float x_former[MAX_n]; int i,j,k; for(i=1;i<=n;++i)
if(fabs(x1[i]-x2[i])>PRECISION) return 1; return 0;
else
}while(IsSatisfyPricision(x,x_former,n) && k