解线性方程组（二）——Jacobi迭代法求解（C++）

迭代法

相比于直接法求解，迭代法使用多次迭代来逐渐逼近解，其精度比不上直接法，但是其速度会比直接法快很多，计算精度可控，特别适用于求解系数矩阵为大型稀疏矩阵的方程组。

Jacobi迭代法

假设有方程组如下：
$$\{\begin{array}{l}{a}_{11}{x}_1+a_{12}{x}_2+\cdots +a_{1n}{x}_n=b_1\\ a_{21}{x}_1+a_{22}{x}_2+\cdots +a_{2n}{x}_n=b_2\\ \cdots \cdots \cdots \\ a_{n1}{x}_1+a_{n2}{x}_2+\cdots +a_{nn}{x}_n=b_n\end{array}$$
将其转换为矩阵形式
$$A\vec{x}=\vec{b}$$
$$\left[\begin{array}{cccc}{a}_{11}& a_{12}& \cdots & a_{1n}\\ a_{21}& a_{22}& \cdots & a_{2n}\\ ⋮& ⋮& \ddots & ⋮\\ a_{m1}& a_{m2}& \cdots & a_{mn}\end{array}\right]\left[\begin{array}{c}{x}_1\\ x_2\\ ⋮\\ x_n\end{array}\right]=\left[\begin{array}{c}{b}_1\\ b_2\\ ⋮\\ b_n\end{array}\right]$$
对于是否可以使用Jacobi迭代法，需要满足以下条件之一：

1. A为行对角优阵，即$|a_{ii}|>{\sum }_{j\ne i}|a_{ij}|(i=1,2,\cdots ,n)$
2. A为行列角优阵，即$|a_{jj}|>{\sum }_{j\ne i}|a_{ij}|(j=1,2,\cdots ,n)$
3. A的元素满足${\sum }_{i\ne j}\frac{|a_{ij}|}{|aii|}<1(j,1,2,\cdots ,n)$
   若矩阵A满足上述条件之一，则可以使用Jacobi迭代法求解方程组。
   首先将上述的方程组转为如下形式：
   $$\{\begin{array}{l}{x}_1=\frac{1}{a_{11}}(-a_{12}{x}_2-\cdots -a_{1n}{x}_n+b_1)\\ x_2=\frac{1}{a_{22}}(-a_{21}{x}_1-\cdots -a_{2n}{x}_n+b_2)\\ \cdots \cdots \cdots \\ x_n=\frac{1}{a_{nn}}(-a_{n1}{x}_1-\cdots -a_{nn-1}{x}_{n-1}+b_n)\end{array}$$
   写成矩阵形式可以得到Jacobi迭代式：
   $$\begin{gathered} (D+L+u)\vec{x}=\vec{b} \\ D\vec{x}=-(L+U)\vec{x}+\vec{b} \\ {\vec{x}}^{(k+1)}=-D^{-1}(L+U){\vec{x}}^{(k)}+D^{-1}\vec{b} \end{gathered}$$
   其中$D$为对角矩阵，$L$为下三角矩阵-$D$，$U$为上三角矩阵-$U$，$D+L+U$为矩阵A。
   

代码实现

由于这个过程涉及大量的矩阵操作，整个算法分为两个源文件：Matrix.cpp实现矩阵操作，main.cpp实现Jacobi迭代法。
首先是Matrix.cpp的代码，其中矩阵求逆的原理参考：

#include <Matrix.h>
#include <iostream>
#include <cmath>
//矩阵与向量相乘，输入矩阵A，向量b，运算结果result和维数n
void matrix_multiply_vector(double **A,double *b,double * result,int n)
{for(int i=0;i<n;i++){result[i]=0.0;for(int j=0;j<n;j++){result[i]+=A[i][j]*b[j];}}
}
//矩阵乘法
void matrix_multiply_matrix(double **A,double **B,double **result,int n)
{for(int i=0;i<n;i++){for(int j=0;j<n;j++){result[i][j]=0.0;for(int k=0;k<n;k++){result[i][j]+=A[i][k]*B[k][j];}}}
}
//矩阵加减法
void matrix_add_matrix(double **A,double **B,double **result,int n,bool isAdd)
{for(int i=0;i<n;i++){for(int j=0;j<n;j++){if(isAdd){result[i][j]=A[i][j]+B[i][j];}else{result[i][j]=A[i][j]-B[i][j];}}}
}
//向量的加减法
void vactor_add_vector(double *A,double *B,double *result,int n,bool isAdd)
{for(int i=0;i<n;i++){if(isAdd){result[i]=A[i]+B[i];}else{result[i]=A[i]-B[i];}}
}
//判断向量误差范围，只要符合精度即可
bool vector_equal(double *A,double *B,int n,double error)
{for(int i=0;i<n;i++){if(fabs(A[i]-B[i])>error){return false;}}return true;
}
//向量赋值
void vector_copy(double *A,double *B,int n)
{for(int i=0;i<n;i++){B[i]=A[i];}
}
//矩阵初始化
void matrix_init(double **A,int n)
{for(int i=0;i<n;i++){A[i]=new double [n];for(int j=0;j<n;j++){A[i][j]=0.0;}}
}
//判断矩阵A是否有收敛性
bool astringency(double **A,int n)
{double abs_row_sum=0.0;double abs_col_sum=0.0;double the_third_condition=0.0;bool RowOptimalMatrix=true;bool ColOptimalMatrix=true;for(int i=0;i<n;i++)//判断是不是行对角优阵{abs_row_sum=0.0;for(int j=0;j<n;j++){if(i!=j){abs_row_sum+=fabs(A[i][j]);}}if(abs_row_sum>A[i][i])//证明不是行对角优阵{RowOptimalMatrix=false;break;}}for(int j=0;j<n;j++)//判断是不是列对角优阵{abs_col_sum=0.0;for(int i=0;i<n;i++){if(i!=j){abs_col_sum+=fabs(A[i][j]);}}if(abs_col_sum>A[j][j]){ColOptimalMatrix=false;break;}}return ColOptimalMatrix or RowOptimalMatrix;
}
//矩阵交换某两行
void matrix_swap_row(double **A,int i,int j,int n)
{double temp;for(int k=0;k<n;k++){temp=A[i][k];A[i][k]=A[j][k];A[j][k]=temp;}
}
//矩阵第i行=矩阵第i行-矩阵第j行*a
void matrix_minus_inner(double **A,double a,int i,int j,int n)
{for(int k=0;k<n;k++){A[i][k]-=a*A[j][k];}
}
//矩阵求逆
void matrix_inverse(double **A,double **A_inverse,int n)
{double **A_E=new double*[2*n];//构建增广矩阵for(int i=0;i<n;i++){A_E[i]=new double [n*2];for(int j=0;j<n*2;j++){if(j<n){A_E[i][j]=A[i][j];}else if((j-n)==i){A_E[i][j]=1;}else{A_E[i][j]=0;}}}//首先将矩阵化为上三角矩阵for(int i=0;i<n;i++){if(A_E[i][i]==0){for(int k=i+1;k<n;k++){if(A_E[k][i]!=0){matrix_swap_row(A_E,i,k,n*2);break;}}}for(int j=i+1;j<n;j++){matrix_minus_inner(A_E,A_E[j][i]/A_E[i][i],j,i,2*n);}}//判断矩阵是否可逆for(int i=0;i<n;i++){if(A_E[i][i]==0){std::cout<<"矩阵不可逆"<<std::endl;exit(0);}}//将上三角转换为对角矩阵for(int j=1;j<n;j++){for(int i=0;i<j;i++){matrix_minus_inner(A_E,A_E[i][j]/A_E[j][j],i,j,2*n);}}for(int i=0;i<n;i++){for(int j=n;j<2*n;j++){A_inverse[i][j-n]=A_E[i][j]/A_E[i][i];}}
}

main.cpp文件内容如下：

//Jacobi迭代法求解线性方程组
/*
5x1+2x2-2x3=1
x1+4x2+x3=2
x1-2x2+4x3=-1
*/
#include<iostream>
#include<cmath>
#include<Matrix.h>//自定义头文件
using namespace std;
int main()
{int n;cout<<"Enter the matrix dimension A: ";cin>>n;//输入数组维度double **A=new double *[n];cout<<"Enter the coefficient matrix:"<<endl;for(int i=0;i<n;i++){A[i]=new double[n];for(int j=0;j<n;j++){cin>>A[i][j];//每次输入一个数字都用空格隔开,输入样例//1 2 3\enter//4 5 6\enter//7 8 9\enter}}double *b=new double[n];cout<<"Input vectors b: ";for(int i=0;i<n;i++){cin>>b[i];//输入方程组右边的向量，1 2 3\enter}bool isAstringency=astringency(A,n);//判断系数矩阵A是否具有收敛性if(isAstringency){cout<<"矩阵A符合收敛性"<<endl;}else{exit(0);cout<<"矩阵A不符合收敛性"<<endl;}double *x=new double[n];//解向量Xdouble *x_last=new double[n];//上一次的xfor(int i=0;i<n;i++){x[i]=0.0;//初始化x}double **A_L_U=new double*[n];//L+Udouble **A_D_inverse=new double*[n];//D的逆for(int i=0;i<n;i++){A_D_inverse[i]=new double [n];A_L_U[i]=new double [n];for(int j=0;j<n;j++){if(i==j){A_L_U[i][j]=0.0;A_D_inverse[i][j]=1.0/A[i][j];//对角矩阵的逆为其倒数}else{A_L_U[i][j]=A[i][j];A_D_inverse[i][j]=0.0;}}}double **B=new double *[n];//公式前半段的矩阵matrix_init(B,n);matrix_multiply_matrix(A_D_inverse,A_L_U,B,n);//求D^(-1)(L+U)double *f=new double[n];matrix_multiply_vector(A_D_inverse,b,f,n);//求取D^-1 * bdouble *temp1=new double[n];do{vector_copy(x,x_last,n);matrix_multiply_vector(B,x_last,temp1,n);//计算公式前半段vactor_add_vector(f,temp1,x,n,false);}while(vector_equal(x,x_last,n,1e-6)==false);//判断向量在误差范围内相等cout<<"运行结果为："<<endl;for(int i=0;i<n;i++){cout<<x[i]<<" ";}system("pause");return 0;
}

结果分析

代码运行结果如下：

当下一次的迭代结果与上一次的迭代结果的最大相差值小于1e-6时，认为迭代已经收敛，输出结果即可（当然也可以换成其它结束迭代方法，如：判断两个向量之差的二范数）。
与直接使用克拉默法则计算准确的解以及matlab计算结果比较，不难发现其$x_1$和$x_3$均不为0，只是是一个在我们设定的误差范围内接近0的数，符合迭代法的解的性质，只能在设定的误差范围内得到一个近似的解。