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

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

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

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

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

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

编译代码：

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

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

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

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

#include
#include
#include
using namespace std;
using namespace Eigen;
const int N=210;
const double eps=1e-6,eeps=1e-6;
int n,flag;
int mp[10][10];
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(auto &acp){		//输出生产的矩阵
	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(char s[10],int tm){
	printf("%s ans is:nx = [ ",s);
	for(int i=1;i<=n;i++) printf("%lf ",x[i][0]);
	printf("]n");
	printf("iteration %d timesn",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(auto &a,auto &b){		//两矩阵相乘
	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(auto &a){		//计算矩阵谱半径
	double mx=0;
	MatrixXd m(n,n);	
	for(int i=0;i es(m,false);
	auto values=es.eigenvalues();
	for(int i=0;ia[mx][i]) mx=j;
		if(abs(a[mx][i])=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 fastern",JNum 
运行结果如下，
    
     
      
       
        ρ
       
       
        (
       
       
        J
       
       
        )
       
       
        =
       
       
        0.84
       
      
      
       rho(J)=0.84
      
     
    ρ(J)=0.84 、
    
     
      
       
        ρ
       
       
        (
       
       
        G
       
       
        )
       
       
        =
       
       
        0.28
       
      
      
       rho(G)=0.28
      
     
    ρ(G)=0.28 。
    
     
      
       
        J
       
       
        a
       
       
        c
       
       
        o
       
       
        b
       
       
        i
       
      
      
       Jacobi
      
     
    Jacobi 法迭代 
    
     
      
       
        80
       
      
      
       80
      
     
    80 次得到答案， 
    
     
      
       
        G
       
       
        a
       
       
        u
       
       
        s
       
       
        s
       
       
        S
       
       
        e
       
       
        i
       
       
        d
       
       
        e
       
       
        l
       
      
      
       GaussSeidel
      
     
    GaussSeidel 法迭代 
    
     
      
       
        15
       
      
      
       15
      
     
    15 次得到答案。
 

  

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

     
      
       
        
         J
        
        
         a
        
        
         c
        
        
         o
        
        
         b
        
        
         i
        
       
       
        Jacobi
       
      
     Jacobi 法与 
     
      
       
        
         G
        
        
         a
        
        
         u
        
        
         s
        
        
         s
        
        
         S
        
        
         e
        
        
         i
        
        
         d
        
        
         e
        
        
         l
        
       
       
        GaussSeidel
       
      
     GaussSeidel 法均不收敛的矩阵：
 

  


     
      
       
        
         J
        
        
         a
        
        
         c
        
        
         o
        
        
         b
        
        
         i
        
       
       
        Jacobi
       
      
     Jacobi 法收敛、 
     
      
       
        
         G
        
        
         a
        
        
         u
        
        
         s
        
        
         s
        
        
         S
        
        
         e
        
        
         i
        
        
         d
        
        
         e
        
        
         l
        
       
       
        GaussSeidel
       
      
     GaussSeidel 法不收敛的矩阵：
 

  


     
      
       
        
         J
        
        
         a
        
        
         c
        
        
         o
        
        
         b
        
        
         i
        
       
       
        Jacobi
       
      
     Jacobi 法不收敛、 
     
      
       
        
         G
        
        
         a
        
        
         u
        
        
         s
        
        
         s
        
        
         S
        
        
         e
        
        
         i
        
        
         d
        
        
         e
        
        
         l
        
       
       
        GaussSeidel
       
      
     GaussSeidel 法收敛的矩阵：
 

  


     
      
       
        
         J
        
        
         a
        
        
         c
        
        
         o
        
        
         b
        
        
         i
        
       
       
        Jacobi
       
      
     Jacobi 法与 
     
      
       
        
         G
        
        
         a
        
        
         u
        
        
         s
        
        
         s
        
        
         S
        
        
         e
        
        
         i
        
        
         d
        
        
         e
        
        
         l
        
       
       
        GaussSeidel
       
      
     GaussSeidel 法均收敛的矩阵：