jacobi2.c

本程序是 用 Jacobi 迭代法解方程组矩阵应

/* 本程序是用Jacobi迭代法解方程组 */
/* 原方程 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]
# define a_jj a[j*n+j]
# define a_ij a[i*n+j]
# define a_ji a[j*n+i]



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 *Init_x(int n)
{
	int i;
	char c;
	double *x;
    x=(double *)malloc(sizeof(double)*n);
	printf("\n是否同意将 X 的初值设为X=(0,0,···,0)\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",&x[i]);
	}
	else for(i=0;i<n;i++) x[i]=0;
	return x;
}

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

void DieDai(double *a,double *b,double *x,int n)
{
    void Out(double *x,int n);
	int i,j;
	long k=0;
    double s,sum,t,jd;
	fflush(stdin);    /* 清除缓冲区的数据 */
	printf("\n请输入迭代终止条件 解的精度 ||Xn-X(n-1)||∞: ");
	scanf("%lf",&jd);
    do
	{
		for(i=0;i<n;i++) a_ii=x[i];
		for(i=0;i<n;i++)
		{
			for(j=0,sum=0;j<n;j++)   /* 求x */
			{
				if(j==i) continue; 
				sum=sum+a_ij*a_jj;
			}
			x[i]=sum+b[i];
		}
		for(i=0,t=0;i<n;i++)    /*  求 ||X||∞   */
		{
			s=doubleAbs(x[i]-a_ii);
			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(x,n);
            printf("★★★请选择 按 0 终止迭代★按 其它任意键 继续迭代下去★★★\n");
			fflush(stdin);     /* 清除缓冲区的数据 */
			if(getch()=='0') break;
		}
	}while(t>jd);
    system("cls");
    printf("\n迭代次数为 %d次",k);
    printf("\n迭代终止时解的精度 ||Xn-X(n-1)||∞=%G\n",t);
}    /*  DieDai end  */

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

void main()
{
	int n;
	char c='1';
	double *a,*b,*x;
   
	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);
		x=Init_x(n);
		DieDai(a,b,x,n);
		Out(x,n);
		free(a);
		free(b);
        free(x);
		printf("\n按 0 退出\n按其它任意键继续解其它方程组"); 
		c=getch();
	}
}