Jacobi迭代法的C++代码实现

Jacobi迭代法，也称为简单迭代法。

假设有方程组

用矩阵表示：

A是系数矩阵，非奇异，b,X为n维向量

假设aii≠0 ，则可以通过移项和简单的乘除法得到如下变换：

其中

构造一个新的矩阵B

则上述经过变换的方程组可以表示为X=BX+g

由此可以构造一个迭代格式（Jacobi迭代）

取定初始向量X(0)之后，便可以由上式产生X(1), X(2), ... X(K), ...

若能够收敛于X(*),则X(*)就是方程的解。

代码实现：

首先定义几个函数用于辅助计算

//矩阵乘法，传入a，b矩阵及其相应的行列数
void mulMatrix(double **&a,int r1,int c1,double **&b,int r2,int c2)
{
	if(c1!=r2) 
		return;
	else { 
		for(int i=0;i<r1;i++)
		{ 
			for(int j=0;j<c1;j++)
			{
				res[i][j]=0; 
				for(int t=0;t<c1;t++)
				res[i][j]=res[i][j]+a[i][t]*b[t][j];}
		}
	}
}

//矩阵加法，传入a,b矩阵及其行列数
void addMatrix(double **&a,double **&b,int r,int c)
{
	for(int i=0;i<r;i++)
		for(int j=0;j<c;j++)
			res[i][j]=a[i][j]+b[i][j];
}

//用于给另外一个矩阵赋值
void moveMatrix(double **&a,double **&b,int r,int c)
{
	for(int i=0;i<r;i++)
		for(int j=0;j<c;j++)
			a[i][j]=b[i][j];
}

//输出矩阵的值
void showMatrix(double **&a,int r,int c)
{
	for(int i=0;i<r;i++)
	{
		for(int j=0;j<c;j++)
			cout<<a[i][j]<<'\t';
		cout<<endl;
	}
}

主函数中并没有去判断结果是否收敛，而是将每一次迭代的结果都打印出来。若要改成判断收敛，只需要判断两次的结果差值是否小于一个可接受的范围。

int main()
{
	int n;
	cout<<"How many rows?"<<endl;
//输入方程变量个数
	cin>>n;
	cout<<"Input matrix:"<<endl;
//输入矩阵A
	double **a=new double*[n];
	for(int i=0;i<n;i++)
		a[i]=new double[n];
	for(int i=0;i<n;i++)
		for(int j=0;j<n;j++)
			cin>>a[i][j];
	cout<<"Calculating matrix B ..."<<endl;
//计算矩阵B
	double **B=new double*[n];
	for(int i=0;i<n;i++)
		B[i]=new double[n];
	//Calculating
	for(int i=0;i<n;i++)
	{
		for(int j=0;j<n;j++)
		{
			if(i!=j)
				B[i][j]=-a[i][j]/a[i][i];
			else 
				B[i][j]=0;
		}
	}
	showMatrix(B,n,n);
	cout<<"Input list b:"<<endl;
//输入方程组右端的n维向量b
	double **b=new double*[n];
	for(int i=0;i<n;i++)
		b[i]=new double[1];
	for(int i=0;i<n;i++)
		for(int j=0;j<1;j++)
			cin>>b[i][j];
	cout<<"Calculating g ..."<<endl;
//计算g
	double **g=new double*[n];
	for(int i=0;i<n;i++)
		g[i]=new double[1];
	//Calculating
	for(int i=0;i<n;i++)
		for(int j=0;j<1;j++)
			g[i][j]=b[i][j]/a[i][i];
	showMatrix(g,n,1);
	//Initial res
//res为全局变量
	res=new double*[n];
	for(int i=0;i<n;i++)
		res[i]=new double[1];
	cout<<"Input X0:"<<endl;
//输入初始向量X(0)
	double **X=new double*[n];
	for(int i=0;i<n;i++)
		X[i]=new double[1];
	for(int i=0;i<n;i++)
		for(int j=0;j<1;j++)
			cin>>X[i][j];
	int k=0;
	cout<<"The result of X"<<k<<":    >>>"<<endl;
	mulMatrix(B,n,n,X,n,1);
	addMatrix(res,g,n,1);
	showMatrix(X,n,1);
//继续计算则按Y，其他则退出计算
	cout<<"Continue?(Y or other else)"<<endl;
	char F;
	cin>>F;
	while(F=='Y'){
		k++;
		cout<<"The result of X"<<k<<":    >>>"<<endl;
		moveMatrix(X,res,n,1);
		mulMatrix(B,n,n,X,n,1);
		addMatrix(res,g,n,1);
		showMatrix(X,n,1);
		cout<<"Continue?(Y/N)"<<endl;
		cin>>F;
	}
}

运行结果：

当计算到X(20)的时候，可以看到已经和X(19)的结果一样了，说明这个方程组用简单迭代法是可以收敛的。

该方程组的准确解为：x1=55/238, x2=5/34, x3=121/238，换算成小数，可以知道上边的计算结果能够收敛到准确解。

算法参考自《数值分析》（华南理工大学出版社），测试样例来自书中的习题。