jacobi.cpp

数值计算方法中几个重要的算法用VC++实现

/*-雅可比(Jacobi)迭代法-*/
#include<stdio.h>
#include<math.h>
#include<stdlib.h>
#define N  10    /*-未知数的最大个数-*/
#define epslion  1e-6
#define maxEps 1e3   /*-根发散的条件-*/
#define type "%lf"
#define Otype "%10.4lf"
typedef double Dtype;

void main()
{
	int i,j,k,n,flag;
	Dtype A[N+1][N+1],b[N+1],x0[N+1],x1[N+1],temp,sum,maxValue0,maxValue1;
	
	printf("Please input n (n<10&&n>0):");
	scanf("%d",&n);
    
	/*-输入系数矩阵A-*/
	printf("Please input Matrix A[%d*%d]:\n",n,n);
    for(i=1;i<=n;i++)
		for(j=1;j<=n;j++)
			scanf(type,&A[i][j]);

	/*-输入向量b-*/
	printf("Please input Vector b[1..%d];\n",n);
	for(i=1;i<=n;i++)
		scanf(type,&b[i]);

	/*-输入初始化向量x(0)-*/
	printf("Please input initialize Vector x0[1..%d]:\n",n);
	for(i=1;i<=n;i++)
		scanf(type,&x0[i]);

	/*-利用雅可比迭代公式计算根x1,x2....xn-*/
	k=0;
	do
	{
		for(i=1;i<=n;i++)
		{
			sum=0;
			for(j=1;j<=n;j++)
				if(j!=i) sum+=A[i][j]*x0[j];
			x1[i]=1/A[i][i]*(b[i]-sum);
		}
		/*-计算||x(k+1)-x(k)||∞-*/
	    if(k>0) maxValue0=maxValue1;
		maxValue1=fabs(x1[1]-x0[1]);
		for(i=2;i<=n;i++)
			if (maxValue1<fabs(x1[i]-x0[i]))
				maxValue1=fabs(x1[i]-x0[i]);
		
		if(maxValue1<epslion) flag=0;
		else 
		{
			flag=1;
			for(i=1;i<=n;i++)
				x0[i]=x1[i];
		}
		if(k>0)
		{
			if(maxValue1-maxValue0>maxEps) /*-发散不收敛，终止程序运行-*/
			{
				printf("Can't use Iterator, The program terminate!");
				exit(0);
			}
		}
		k=k+1;
	}while(flag);
	/*-输出x1,x2,x3....xn-*/
	printf("The result is :\n");
	for(i=1;i<=n;i++)
		printf("x[%d]="Otype"\n",i,x1[i]);

}