迭代法
相比于直接法求解,迭代法使用多次迭代来逐渐逼近解,其精度比不上直接法,但是其速度会比直接法快很多,计算精度可控,特别适用于求解系数矩阵为大型稀疏矩阵的方程组。
Jacobi迭代法
假设有方程组如下:
{ a 11 x 1 + a 12 x 2 + ⋯ + a 1 n x n = b 1 a 21 x 1 + a 22 x 2 + ⋯ + a 2 n x n = b 2 ⋯ ⋯ ⋯ a n 1 x 1 + a n 2 x 2 + ⋯ + a n n x n = b n \begin{cases} 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 \qquad \qquad\cdots \qquad \qquad \cdots \\ a_{n1}x_1+a_{n2}x_2+\cdots+a_{nn}x_n=b_n\\ \end{cases} ⎩ ⎨ ⎧a11x1+a12x2+⋯+a1nxn=b1a21x1+a22x2+⋯+a2nxn=b2⋯⋯⋯an1x1+an2x2+⋯+annxn=bn
将其转换为矩阵形式
A x ⃗ = b ⃗ A\vec{x}=\vec{b} Ax=b
[ a 11 a 12 ⋯ a 1 n a 21 a 22 ⋯ a 2 n ⋮ ⋮ ⋱ ⋮ a m 1 a m 2 ⋯ a m n ] [ x 1 x 2 ⋮ x n ] = [ b 1 b 2 ⋮ b n ] \begin{bmatrix} {a_{11}}&{a_{12}}&{\cdots}&{a_{1n}}\\ {a_{21}}&{a_{22}}&{\cdots}&{a_{2n}}\\ {\vdots}&{\vdots}&{\ddots}&{\vdots}\\ {a_{m1}}&{a_{m2}}&{\cdots}&{a_{mn}}\\ \end{bmatrix} \begin{bmatrix} {x_{1}}\\ {x_{2}}\\ {\vdots}\\ {x_{n}}\\ \end{bmatrix}= \begin{bmatrix} {b_{1}}\\ {b_{2}}\\ {\vdots}\\ {b_n} \end{bmatrix} a11a21⋮am1a12a22⋮am2⋯⋯⋱⋯a1na2n⋮amn x1x2⋮xn = b1b2⋮bn
对于是否可以使用Jacobi迭代法,需要满足以下条件之一:
- A为行对角优阵,即 ∣ a i i ∣ > ∑ j ≠ i ∣ a i j ∣ ( i = 1 , 2 , ⋯ , n ) |a_{ii}|>\sum_{j \neq i}|a_{ij}|(i=1,2,\cdots,n) ∣aii∣>∑j=i∣aij∣(i=1,2,⋯,n)
- A为行列角优阵,即 ∣ a j j ∣ > ∑ j ≠ i ∣ a i j ∣ ( j = 1 , 2 , ⋯ , n ) |a_{jj}|>\sum_{j \neq i}|a_{ij}|(j=1,2,\cdots,n) ∣ajj∣>∑j=i∣aij∣(j=1,2,⋯,n)
- A的元素满足 ∑ i ≠ j ∣ a i j ∣ ∣ a i i ∣ < 1 ( j , 1 , 2 , ⋯ , n ) \sum_{i \neq j}\frac{|a_{ij}|}{|aii|}<1(j,1,2,\cdots,n) ∑i=j∣aii∣∣aij∣<1(j,1,2,⋯,n)
若矩阵A满足上述条件之一,则可以使用Jacobi迭代法求解方程组。
首先将上述的方程组转为如下形式:
{ x 1 = 1 a 11 ( − a 12 x 2 − ⋯ − a 1 n x n + b 1 ) x 2 = 1 a 22 ( − a 21 x 1 − ⋯ − a 2 n x n + b 2 ) ⋯ ⋯ ⋯ x n = 1 a n n ( − a n 1 x 1 − ⋯ − a n n − 1 x n − 1 + b n ) \begin{cases} 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 \qquad \qquad\cdots \qquad \qquad \cdots \\ x_n=\frac{1}{a_{nn}}(-a_{n1}x_1-\cdots -a_{nn-1}x_{n-1}+b_n)\\ \end{cases} ⎩ ⎨ ⎧x1=a111(−a12x2−⋯−a1nxn+b1)x2=a221(−a21x1−⋯−a2nxn+b2)⋯⋯⋯xn=ann1(−an1x1−⋯−ann−1xn−1+bn)
写成矩阵形式可以得到Jacobi迭代式:
( D + L + u ) x ⃗ = b ⃗ D x ⃗ = − ( L + U ) x ⃗ + b ⃗ x ⃗ ( k + 1 ) = − D − 1 ( L + U ) x ⃗ ( k ) + D − 1 b ⃗ (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} (D+L+u)x=bDx=−(L+U)x+bx(k+1)=−D−1(L+U)x(k)+D−1b
其中 D D D为对角矩阵, L L L为下三角矩阵- D D D, U U U为上三角矩阵- U U U, D + L + U D+L+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_1 x1和 x 3 x_3 x3均不为0,只是是一个在我们设定的误差范围内接近0的数,符合迭代法的解的性质,只能在设定的误差范围内得到一个近似的解。