雅可比迭代与高斯-赛德尔迭代的C++实现，及判断收敛性

随机生成若干个 $n$ 阶方阵与 $n$ 阶向量构成 $Ax=b$

• 分别判断J法和GS法的收敛性
  • 是否能收敛
    • 报告中应生成部分不收敛的矩阵
  • 估计J法和GS法收敛速度哪个更快
• 实现用J法和GS法解该方程组
• 实验判断J法和GS法的收敛速度，并与理论估计作对比

首先随机生成矩阵，判断它是否非奇异、判断对角矩阵D是否非奇异、判断J法和GS法是否都收敛，如果不满足就重新生成。复杂度很高

得到矩阵后，通过 $\rho (J)$ 与 $\rho (G)$ 的大小可以估计哪种迭代法收敛速度更快。然后使用 J 法和 GS 法分别求解，当误差小于给定值时，退出循环。记录迭代次数，就可以验证之前的估计是否正确。

有了上次实验中计算矩阵最大特征值的经历后，这次计算矩阵谱半径就容易很多了。注意实矩阵的特征值可能为复数，复数的模等于复平面上向量的模。不要只计算实部，否则生成的矩阵有概率不收敛，求出来的值都是 nan 与 inf 。

编译代码：

g++ main.cpp -o main -I .\eigen3 -w -O2

其中 .\eigen3 为第三方库 Eigen 的相对路径， -w 忽略警告，-O2 开启编译优化。
Eigen 库的安装使用教程： 链接

程序的功能为：输入矩阵阶数，随机生成一个 J 法与 GS 法均收敛的矩阵，求出 $\rho (J)$ 与 $\rho (G)$ ，并且以此估计两种方法的收敛速度。然后分别用这两种方法求解，并输出达到预定精度所需的迭代次数。

命令行运行，否则程序运行结束后，控制台不会停留。

#include<bits/stdc++.h>
#include<Eigen/Dense>
#include<Eigen/Eigenvalues>
using namespace std;
using namespace Eigen;
const int N=210;
const double eps=1e-6,eeps=1e-6;
int n,flag;
double a[N][N],acp[N][N],tmp[N][N],t[N][N],x[N][2];

void randInit(){	//随机生成矩阵
	for(int i=1;i<=n;i++) for(int j=1;j<=n+1;j++)
		acp[i][j]=rand()%200/10.0-10;
	for(int i=1;i<=n;i++) if(acp[i][i]==0)
		acp[i][i]=(rand()%100/10.0+1)*((rand()&1)?1:-1);
}

void out(double acp[N][N]){		//输出生产的矩阵
	puts("A is:");
	for(int i=1;i<=n;i++) for(int j=1;j<=n+1;j++)
		printf("%lf%c",acp[i][j]," \n"[j==n+1]);
}

void pt(string s,int tm){
	printf("%s ans is:\nx = [ ",s.c_str());
	for(int i=1;i<=n;i++) printf("%lf ",x[i][0]);
	printf("]\n");
	printf("iteration %d times\n",tm);
}

void init(){		//初始化，把生成的矩阵复制到A
	for(int i=1;i<=n;i++) for(int j=1;j<=n+1;j++)
		a[i][j]=acp[i][j];
}

void mul(double a[N][N],double b[N][N]){		//两矩阵相乘
	for(int i=1;i<=n;i++)
		for(int j=1;j<=n;j++){
			tmp[i][j]=0;
			for(int k=1;k<=n;k++)
				tmp[i][j]+=a[i][k]*b[k][j];
		}
	for(int i=1;i<=n;i++)
		for(int j=1;j<=n;j++)
			a[i][j]=tmp[i][j];
}

double calc(double a[N][N]){		//计算矩阵谱半径
	double mx=0;
	MatrixXd m(n,n);	
	for(int i=0;i<n;i++)
		for(int j=0;j<n;j++)
			m(i,j)=a[i+1][j+1];
	EigenSolver<MatrixXd> es(m,false);
	auto values=es.eigenvalues();
	for(int i=0;i<n;i++)
		mx=max(mx,sqrt(values[i].real()*values[i].real()+values[i].imag()*values[i].imag()));
	return mx;
}

void getReverse(double a[N][N]){	//矩阵求逆
	for(int i=1;i<=n;i++) for(int j=n+1;j<=n+n;j++) a[i][j]=(i+n==j);
	for(int i=1;i<=n;i++){
		int mx=i;
		for(int j=i+1;j<=n;j++)
			if(a[j][i]>a[mx][i]) mx=j;
		if(abs(a[mx][i])<eps) return (void)(flag=true);
		for(int j=i;j<=n+n;j++) swap(a[i][j],a[mx][j]);
		for(int j=n+n;j>=i;j--) a[i][j]/=a[i][i];
		for(int k=1;k<=n;k++) if(i!=k)
			for(int j=n+n;j>=i;j--) a[k][j]-=a[k][i]*a[i][j];
	}
	for(int i=1;i<=n;i++) for(int j=1;j<=n;j++) a[i][j]=a[i][j+n];
}

double judgeJ(){	//判断J法是否收敛
	init();	
	for(int i=1;i<=n;i++) for(int j=1;j<=n;j++){
		t[i][j]=0;
		if(i==j){
			t[i][j]=1.0/a[i][j];
			a[i][j]=0;
		}
		a[i][j]=-a[i][j];
	}
	mul(t,a);
	return calc(t);
}

double judgeGS(){		//判断GS法是否收敛
	init();
	for(int i=1;i<=n;i++) for(int j=1;j<=n;j++)
		if(j>i) t[i][j]=-a[i][j],a[i][j]=0;
		else t[i][j]=0;
	getReverse(a);
	mul(a,t);
	return calc(a);
}

int Jacobi(){	//雅可比迭代法
	int cur=0,tm=0;
	double mx;
	for(int i=1;i<=n;i++) x[i][0]=0;
	do{
		tm++; mx=0;
		for(int i=1;i<=n;i++){
			x[i][cur^1]=a[i][n+1];
			for(int j=1;j<=n;j++){
				if(j==i) continue;
				x[i][cur^1]-=a[i][j]*x[j][cur];
			}
			x[i][cur^1]/=a[i][i];
		} 
		for(int i=1;i<=n;i++) mx=max(mx,abs(x[i][0]-x[i][1]));
		cur^=1;
	}while(mx>eeps);
	for(int i=1;i<=n;i++) x[i][0]=x[i][cur];
	return tm;
}

int GaussSeidel(){	//高斯-赛德尔迭代法
	init();
	for(int i=1;i<=n;i++) x[i][0]=0;
	int tm=0;
	double mx;
	do{
		tm++; mx=0;
		for(int i=1;i<=n;i++){
			double t=x[i][0];
			x[i][0]=a[i][n+1];
			for(int j=1;j<=n;j++){
				if(i==j) continue;
				x[i][0]-=a[i][j]*x[j][0];
			}
			x[i][0]/=a[i][i];
			mx=max(mx,abs(x[i][0]-t));
		}
	}while(mx>eeps);
	return tm;
}

int main(){
	double JNum,GSNum;
	srand(time(0));
	puts("input n (2<=n<=10):");
	scanf("%d",&n);
	do{
		randInit();
		flag=false;
		init();
		getReverse(a);	//判断矩阵是否奇异
	}while(!((JNum=judgeJ())<1&&(GSNum=judgeGS())<1)||flag);
	out(acp);
	printf("JNum is %lf, and GSNum is %lf.\n",JNum,GSNum);
	printf("So maybe %s is faster\n",JNum<GSNum?"J":"GS");
	init(); pt("Jacobi",Jacobi());
	init(); pt("GaussSeidel",GaussSeidel());
}

运行结果如下，$\rho (J)=0.84$ 、$\rho (G)=0.28$ 。$Jacobi$ 法迭代 $80$ 次得到答案， $GaussSeidel$ 法迭代 $15$ 次得到答案。

最后生成一些不收敛的矩阵：

1. $Jacobi$ 法与 $GaussSeidel$ 法均不收敛的矩阵：

2. $Jacobi$ 法收敛、 $GaussSeidel$ 法不收敛的矩阵：

3. $Jacobi$ 法不收敛、 $GaussSeidel$ 法收敛的矩阵：

4. $Jacobi$ 法与 $GaussSeidel$ 法均收敛的矩阵：