gauss_seidel.c

guass_sedal法解方程

 
/* 本程序是用 Gauss-Seidel 迭代法解方程组 */
/* 原方程 AX=b 迭代方程 X(k+1)=BoX(k)+f */
/* Bo=I-(D-1)A ；f=(D-1)b. (其中（D-1）表示D的逆） */

# include <process.h> 
# include <stdio.h>
# include <stdlib.h>
# include <conio.h>
# define a_ii a[i*n+i]  /* a_ii就是a[i][i]   */
# define a_jj a[j*n+j]   /* a_jj就是a[j][j]  */
# define a_ij a[i*n+j]   /* a_ij就是a[i][j]  */
# define a_ji a[j*n+i]   /* a_ji就是a[j][i]  */

void  ShuoMing()
{
  system("cls");
  printf("\n\t* * * * * * * * * * * * * * * * * * * * * * *\n");
  printf("\t*                                           *\n");
  printf("\t*                 使用说明                  *\n");
  printf("\t*  =======================================  *\n");
  printf("\t*                                           *\n");
  printf("\t* 本程序是                                  *\n");
  printf("\t*        用 Gauss-Seide 迭代法解方程组      *\n");
  printf("\t*   矩阵应是满足Gauss-Seide迭代法收敛       *\n");
  printf("\t*                                           *\n");
  printf("\t* 使用时： 系数矩阵 与 常数项 分开输入      *\n");
  printf("\t*   输入系数矩阵时只能按行的顺序输入        *\n");
  printf("\t*   即依次按a11,a12...a1n顺序输完第一行     *\n");
  printf("\t*   再依次输入第二行 依次类推直到第n行      *\n");
  printf("\t*   提示输入迭代终止条件时输入的数据        *\n");
  printf("\t*   必需大约或等于0                         *\n");
  printf("\t* 更多详情参见 使用说明文档                 *\n");
  printf("\t*                                           *\n");
  printf("\t* 作者 E-mail [e-mail]   QQ 649990777 *\n");
  printf("\t* * * * * * * * * * * * * * * * * * * * * * *\n");
  printf("\t\t请按任意键继续:");                 
  getch();system("cls");
}

double *Get_A(int n)
{
 	int i,j;
    double *a;
	a=(double *)malloc(sizeof(double)*n*n);
    if(a==NULL)
	{
		printf("给系数矩阵分配存储空间失败 按任意建结束程序\n"); 
		getch();
		exit(0);
	}  
    printf("\n请输入系数矩阵（按行的顺序）\n");
	fflush(stdin);     /* 清除缓冲区的数据 */
	for(i=0;i<n;i++)
		for(j=0;j<n;j++)
			scanf("%lf",&a_ij); 
	return a;
}       /*   Get_A End   */

double *Get_b(int n)
{
	int i;
    double *b;
    b=(double *)malloc(sizeof(double)*n);
	if(b==NULL)
	{
		printf("给常数项分配存储空间失败 按任意建结束程序\n"); 
		getch();
		exit(0);
	}
	printf("请依次输入常数项\n");
	fflush(stdin);    /* 清除缓冲区的数据 */
	for(i=0;i<n;i++)
		scanf("%lf",&b[i]); 
	return b;
}/*   Get_b End   */

void Qiou_f_Bo(double *a,double *b,int n)
{
	int i,j;
    for(i=0;i<n;i++)    /*  求f */
		b[i]=b[i]/a_ii; 
	for(i=0;i<n;i++)
	{
		for(j=0;j<n;j++)       /*  求Bo */
		{
			if(j==i) continue; 
			a_ij=-a_ij/a_ii; 
		}
	}
}     /*  Qiou_f_Bo End */

double doubleAbs(double y)
{
	if(y<0) return -y;
	else return y;
}

void GaussSeidel_DieDai(double *a,double *b,int n)
{
    void Out(double *a,int n);
	char c;
	int i,j;
	long k=0;
    double s,sum,t,jd;
	printf("\n是否同意将 X 的初值设为X=(0,0,···,0)T\n");
    printf(" 不同意：按 1 重新手动输入赋初值(需输入 %d 个值)\n",n);
    printf(" 同意：  按其它任意键确认\n");
	c=getch();
	if(c=='1')
	{
		printf(" 请按x1到xn的顺序输入xi的值\n");
		fflush(stdin);    /* 清除缓冲区的数据 */
		for(i=0;i<n;i++)
			scanf("%lf",&a_ii);
	}
	else for(i=0;i<n;i++) a_ii=0;
	fflush(stdin);         /* 清除缓冲区的数据 */
	printf("\n请输入迭代终止条件 解的精度 ||Xn-X(n-1)||∞: ");
	scanf("%lf",&jd);
    do                           /*  迭代 */
	{
		for(i=0,t=0;i<n;i++)
		{
			for(j=0,sum=0;j<n;j++)   /* 求x */
			{
				if(j==i) continue; 
				sum=sum+a_ij*a_jj;
			}
			s=a_ii;
			a_ii=sum+b[i];
			s=doubleAbs(a_ii-s);  /*  求 ||X||∞  */
			t=s>t?s:t;
		} 
        if(++k%10000000==0)   
		{
			printf("\n\n以经迭代了 %ld 次为防止死循环迭代被强行暂停\n",k);
            printf("当前解的精度 未达到要求 情况如下\n");
            printf("★★★★★当前 解的精度 ||Xn-X(n-1)||∞=%G\n",t);
            Out(a,n);
            printf("★★★请选择 按 0 终止迭代★按 其它任意键 继续迭代下去★★★\n");
			fflush(stdin);         /* 清除缓冲区的数据 */
			if(getch()=='0') break;
		}
	}while(t>jd);
    printf("\n迭代次数为 %d次",k);
    printf("\n迭代终止时解的精度 ||Xn-X(n-1)||∞=%G\n",t);
}    /*  DieDai end  */

void Out(double *a,int n)
{
	int i;
    printf("方程的解为:\n");
	for(i=0;i<n;i++)
		printf("x[%d]=%G \n",i+1,a_ii);
}

void main()
{
	int n;
	char c='1';
	double *a,*b;
    ShuoMing();
	while(c!='0')
	{
		system("cls");
		printf("\n请输入矩阵的阶数(若阶数为 0 则退出）: ");  
		scanf("%d",&n);
		if(n<=0) exit(0);
		a=Get_A(n);     /* 构建并载入矩阵 */
		b=Get_b(n);    /* 构建并载入常数项 */
		Qiou_f_Bo(a,b,n);
		GaussSeidel_DieDai(a,b,n);
		Out(a,n);
		free(a);
		free(b);
		printf("\n按 0 退出\n按 其它任意键 继续解其它方程组"); 
		c=getch();
	}
}