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

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

建议直接生成对称正定矩阵

取 0<ω<20<\omega<20<ω<2 的若干值进行收敛代数对比

用充分条件判断该矩阵 AAA 是否对 SOR 迭代法收敛

实现用 SOR 法解该方程组
与精确解对比并计算误差

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

对称阵的三角分解定理：设 AAA 为 nnn 阶对称矩阵，且 AAA 的所有顺序主子式均不为零，则 AAA 可以唯一分解为 A=LDLTA=LDL^TA=LDLT.
其中 LLL 为单位下三角矩阵，DDD 为对角矩阵，且 DDD 的元素 did_idi 均为正数。

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

放一张程序运行截图：

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

代码内部迭代终止条件为 max1≤i≤n∣xi(k+1)−xi(k)∣<εmax_{1 \le i \le n}|x_i^{(k+1)}-x_i^{(k)}|<\varepsilonmax1≤i≤n∣xi(k+1)−xi(k)∣<ε .

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

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

#include<bits/stdc++.h>
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])<eps) return (void)(flag=true);for(int j=i;j<=n+1;j++) swap(a[i][j],a[mx][j]);for(int j=n+1;j>=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++){       //枚举 ww=i/10.0;printf("w is %.2lf, iterator times is %d, absolute error is %.10lf.\n",w,sor(),absError);out();           //输出解}
}

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

逐次超松弛迭代法SOR
%SOR:逐次超松弛迭代法 %本函数只能求解当A为n*n的矩阵 %2010-10-21 function x=SOR(A,b,w) tic; if nargin==2 ...
sor松弛迭代matlab,数值分析Python实现系列—— 二、逐次超松弛迭代法(SOR)
二.超松弛迭代法(SOR) 1.原理: 回顾: 在一般情况下 : 收敛过慢甚至不收敛的$B$与$f$,经过对系数矩阵$A$分裂成$A = M - N$的形式, 使得迭代公式变为: $x^{k+ ...
3.3 数值分析: 逐次超松弛迭代法-SOR方法
本文内容为东北大学数值分析国家精品慕课课程的课程讲义, 将其整理为OneNote笔记同时添加了本人上课时的课堂笔记, 且主页中的思维导图就是根据课件内容整理而来, 为了方便大家和自己查看,特将此上传到 ...
松弛迭代法matlab,逐次超松弛迭代法求解线性方程组的MATLAB实现
function [X_reality,n_reality] = SOR(A,b,X_start,w,n_limit,tolerance) %% % A为迭代的系数矩阵 % b为方程组右边的常数项(列 ...
MATLAB 五对角矩阵 Jacobi迭代法 SOR迭代法解方程组
% ----------------------------------------------------------------------------------------------- % ...
三种迭代法解方程组（雅可比Jacobi、高斯-赛德尔Gaisi_saideer、逐次超松弛SOR）
分析用下列迭代法解线性方程组 4 -1 0 -1 0 0 0 -1 4 -1 0 -1 0 5 0 -1 4 -1 0 -1 -2 -1 0 -1 4 -1 0 ...
高斯赛德尔迭代c语言_逐次超松弛SOR迭代概述
本文使用 Zhihu On VSCode 创作并发布逐次超松弛(SOR)迭代法概述一.方法背景逐次超松弛迭代法是高斯-塞德尔迭代法的一种变种,是为了解决线性方程组的一种迭代方法.由David ...
超松弛迭代法（SOR）的Python实现
在Gauss-Seidel迭代法中M=D-L的基础上,我们令M=(1/w)(D-wL),B=I-M的逆A,f=M的逆b,其中w为可选择的松弛因子,要求w>0 于是M的逆=w((D-wL)的逆), ...
线性方程组的SOR迭代法
写在前面 1 SOR简介 2 SOR算法推导 3 SOR算法收敛性 4 实例分析 5 代码实现 6 参考文献与链接迭代法收敛定理定理1(充要条件): 矩阵 M \mathbf{M} M的谱半径 ρ ...

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

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

最新文章

热门文章