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*
;
for(int i=0;i<n;i++)
a[i]=new double
;
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*
;
for(int i=0;i<n;i++)
B[i]=new double
;
//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*
;
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*
;
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*
;
for(int i=0;i<n;i++)
res[i]=new double[1];
cout<<"Input X0:"<<endl;
//输入初始向量X(0)
double **X=new double*
;
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，换算成小数，可以知道上边的计算结果能够收敛到准确解。

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