逐次超松弛迭代法 ( SOR ) 的C++实现

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

• 建议直接生成对称正定矩阵
• 取 0 < ω < 2 0
• 用充分条件判断该矩阵 A A A 是否对 SOR 迭代法收敛
• 实现用 SOR 法解该方程组
  • 与精确解对比并计算误差

首先生成一个对称正定矩阵。由对称阵的三角分解定理，先随机生成一个下三角矩阵 L，再生成一个元素均为正数的对角矩阵 D D D，计算 L D L T LDL^T LDLT 就可以得到一个对称正定矩阵。

对称阵的三角分解定理：设 A A A 为 n n n 阶对称矩阵，且 A A A 的所有顺序主子式均不为零，则 A A A 可以唯一分解为 A = L D L T A=LDL^T A=LDLT.
其中 L L L 为单位下三角矩阵， D D D 为对角矩阵，且 D D D 的元素 d i d_i di​ 均为正数。

通过高斯消元，先算一个精确答案出来。再枚举 ω omega ω 的值，调用 SOR 迭代法。SOR 迭代法内部统计迭代次数与绝对误差，将这些信息都输出出来。

放一张程序运行截图：

可以看到，在这个例子中，当 ω omega ω 为 1.6 1.6 1.6 的时候，迭代次数最少。

代码内部迭代终止条件为 m a x 1 ≤ i ≤ n ∣ x i ( k + 1 ) − x i ( k ) ∣ < ε max_{1 le i le n}|x_i^{(k+1)}-x_i^{(k)}|

ε varepsilon ε 的取值会影响迭代次数，以及误差。

程序刚运行时，需要输入矩阵的阶数。不要输入太大的数，否则要运行很久。

#include
using namespace std;
const int N=210;
const double eeps=1e-6,eps=-1e-6;
int n,flag;
double acp[N][N],a[N][N],l[N][N],tmp[N][N],x[N],w,stdX[N],absError;

void init(){ for(int i=1;i<=n;i++) for(int j=1;j<=n+1;j++) a[i][j]=acp[i][j]; }
void out(){ printf("x = [ "); for(int i=1;i<=n;i++) printf("%lf ",x[i]); puts("]"); }
void pt(){ 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 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];
}

void randinit(){	//随机生成对称正定矩阵
	for(int i=1;i<=n;i++)
		for(int j=1;j<=n;j++)
			if(j>i) acp[i][j]=l[j][i]=0;
			else acp[i][j]=l[j][i]=(rand()%100+1)/10.0*(rand()&1?-1:1);
	for(int i=1;i<=n;i++)
		for(int j=1;j<=n;j++)
			if(i==j) a[i][j]=(rand()%10+1);
	mul(acp,a);
	mul(acp,l);
	for(int i=1;i<=n;i++) acp[i][n+1]=(rand()%200)/10.0-10;
}

void gauss(){		//高斯消元计算标准答案
	init();
	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])=i;j--) a[i][j]/=a[i][i];
		for(int k=1;k<=n;k++) if(i!=k)
			for(int j=n+1;j>=i;j--) a[k][j]-=a[k][i]*a[i][j];
	}
	for(int i=1;i<=n;i++) stdX[i]=a[i][n+1];
}

int sor(){		//sor 迭代
	init();
	for(int i=1;i<=n;i++) x[i]=0;		
	double mx;
	int tm=0;
	do{
		tm++; mx=0;
		for(int i=1;i<=n;i++){
			double add=a[i][n+1];
			for(int j=1;j<=n;j++) add-=a[i][j]*x[j];
			add*=w/a[i][i];
			mx=max(mx,abs(add));
			x[i]+=add;
		} 
	}while(mx>eeps);	//迭代结束条件
	absError=0;			//统计绝对误差
	for(int i=1;i<=n;i++) absError+=abs(x[i]-stdX[i]);
	return tm;			//返回迭代次数
}

int main(){
	srand(time(0));
	puts("input n (2<=n<=20):");
	scanf("%d",&n);
	do{
		flag=false;
		randinit();		//随机生成对称正定矩阵
		gauss();		//判断是否非奇异，得出标准答案
	}while(flag);
	pt();				//输出生成的矩阵
	for(int i=1;i<20;i++){		//枚举 w
		w=i/10.0;
		printf("w is %.2lf, iterator times is %d, absolute error is %.10lf.n",w,sor(),absError);
		out();			//输出解
	}
}