免费元素素材网站,安徽网络建站,php模板网站,北京工程造价信息网迭代法
相比于直接法求解#xff0c;迭代法使用多次迭代来逐渐逼近解#xff0c;其精度比不上直接法#xff0c;但是其速度会比直接法快很多#xff0c;计算精度可控#xff0c;特别适用于求解系数矩阵为大型稀疏矩阵的方程组。
Jacobi迭代法
假设有方程组如下#xf…迭代法
相比于直接法求解迭代法使用多次迭代来逐渐逼近解其精度比不上直接法但是其速度会比直接法快很多计算精度可控特别适用于求解系数矩阵为大型稀疏矩阵的方程组。
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_1a_{12}x_2\cdotsa_{1n}x_nb_1\\ a_{21}x_1a_{22}x_2\cdotsa_{2n}x_nb_2\\ \cdots \qquad \qquad\cdots \qquad \qquad \cdots \\ a_{n1}x_1a_{n2}x_2\cdotsa_{nn}x_nb_n\\ \end{cases} ⎩ ⎨ ⎧a11x1a12x2⋯a1nxnb1a21x1a22x2⋯a2nxnb2⋯⋯⋯an1x1an2x2⋯annxnbn 将其转换为矩阵形式 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}|(i1,2,\cdots,n) ∣aii∣∑ji∣aij∣(i1,2,⋯,n)A为行列角优阵即 ∣ a j j ∣ ∑ j ≠ i ∣ a i j ∣ ( j 1 , 2 , ⋯ , n ) |a_{jj}|\sum_{j \neq i}|a_{ij}|(j1,2,\cdots,n) ∣ajj∣∑ji∣aij∣(j1,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) ∑ij∣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_nb_1)\\ x_2\frac{1}{a_{22}}(-a_{21}x_1-\cdots -a_{2n}x_nb_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} ⎩ ⎨ ⎧x1a111(−a12x2−⋯−a1nxnb1)x2a221(−a21x1−⋯−a2nxnb2)⋯⋯⋯xnann1(−an1x1−⋯−ann−1xn−1bn) 写成矩阵形式可以得到Jacobi迭代式 ( D L u ) x ⃗ b ⃗ D x ⃗ − ( L U ) x ⃗ b ⃗ x ⃗ ( k 1 ) − D − 1 ( L U ) x ⃗ ( k ) D − 1 b ⃗ (DLu)\vec{x}\vec{b}\\ D\vec{x}-(LU)\vec{x}\vec{b}\\ \vec{x}^{(k1)}-D^{-1}(LU)\vec{x}^{(k)}D^{-1}\vec{b} (DLu)x b Dx −(LU)x b x (k1)−D−1(LU)x (k)D−1b 其中 D D D为对角矩阵 L L L为下三角矩阵- D D D U U U为上三角矩阵- U U U D L U DLU DLU为矩阵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 i0;in;i){result[i]0.0;for(int j0;jn;j){result[i]A[i][j]*b[j];}}
}
//矩阵乘法
void matrix_multiply_matrix(double **A,double **B,double **result,int n)
{for(int i0;in;i){for(int j0;jn;j){result[i][j]0.0;for(int k0;kn;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 i0;in;i){for(int j0;jn;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 i0;in;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 i0;in;i){if(fabs(A[i]-B[i])error){return false;}}return true;
}
//向量赋值
void vector_copy(double *A,double *B,int n)
{for(int i0;in;i){B[i]A[i];}
}
//矩阵初始化
void matrix_init(double **A,int n)
{for(int i0;in;i){A[i]new double [n];for(int j0;jn;j){A[i][j]0.0;}}
}
//判断矩阵A是否有收敛性
bool astringency(double **A,int n)
{double abs_row_sum0.0;double abs_col_sum0.0;double the_third_condition0.0;bool RowOptimalMatrixtrue;bool ColOptimalMatrixtrue;for(int i0;in;i)//判断是不是行对角优阵{abs_row_sum0.0;for(int j0;jn;j){if(i!j){abs_row_sumfabs(A[i][j]);}}if(abs_row_sumA[i][i])//证明不是行对角优阵{RowOptimalMatrixfalse;break;}}for(int j0;jn;j)//判断是不是列对角优阵{abs_col_sum0.0;for(int i0;in;i){if(i!j){abs_col_sumfabs(A[i][j]);}}if(abs_col_sumA[j][j]){ColOptimalMatrixfalse;break;}}return ColOptimalMatrix or RowOptimalMatrix;
}
//矩阵交换某两行
void matrix_swap_row(double **A,int i,int j,int n)
{double temp;for(int k0;kn;k){tempA[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 k0;kn;k){A[i][k]-a*A[j][k];}
}
//矩阵求逆
void matrix_inverse(double **A,double **A_inverse,int n)
{double **A_Enew double*[2*n];//构建增广矩阵for(int i0;in;i){A_E[i]new double [n*2];for(int j0;jn*2;j){if(jn){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 i0;in;i){if(A_E[i][i]0){for(int ki1;kn;k){if(A_E[k][i]!0){matrix_swap_row(A_E,i,k,n*2);break;}}}for(int ji1;jn;j){matrix_minus_inner(A_E,A_E[j][i]/A_E[i][i],j,i,2*n);}}//判断矩阵是否可逆for(int i0;in;i){if(A_E[i][i]0){std::cout矩阵不可逆std::endl;exit(0);}}//将上三角转换为对角矩阵for(int j1;jn;j){for(int i0;ij;i){matrix_minus_inner(A_E,A_E[i][j]/A_E[j][j],i,j,2*n);}}for(int i0;in;i){for(int jn;j2*n;j){A_inverse[i][j-n]A_E[i][j]/A_E[i][i];}}
}main.cpp文件内容如下
//Jacobi迭代法求解线性方程组
/*
5x12x2-2x31
x14x2x32
x1-2x24x3-1
*/
#includeiostream
#includecmath
#includeMatrix.h//自定义头文件
using namespace std;
int main()
{int n;coutEnter the matrix dimension A: ;cinn;//输入数组维度double **Anew double *[n];coutEnter the coefficient matrix:endl;for(int i0;in;i){A[i]new double[n];for(int j0;jn;j){cinA[i][j];//每次输入一个数字都用空格隔开,输入样例//1 2 3\enter//4 5 6\enter//7 8 9\enter}}double *bnew double[n];coutInput vectors b: ;for(int i0;in;i){cinb[i];//输入方程组右边的向量1 2 3\enter}bool isAstringencyastringency(A,n);//判断系数矩阵A是否具有收敛性if(isAstringency){cout矩阵A符合收敛性endl;}else{exit(0);cout矩阵A不符合收敛性endl;}double *xnew double[n];//解向量Xdouble *x_lastnew double[n];//上一次的xfor(int i0;in;i){x[i]0.0;//初始化x}double **A_L_Unew double*[n];//LUdouble **A_D_inversenew double*[n];//D的逆for(int i0;in;i){A_D_inverse[i]new double [n];A_L_U[i]new double [n];for(int j0;jn;j){if(ij){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 **Bnew double *[n];//公式前半段的矩阵matrix_init(B,n);matrix_multiply_matrix(A_D_inverse,A_L_U,B,n);//求D^(-1)(LU)double *fnew double[n];matrix_multiply_vector(A_D_inverse,b,f,n);//求取D^-1 * bdouble *temp1new 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 i0;in;i){coutx[i] ;}system(pause);return 0;
}结果分析
代码运行结果如下
当下一次的迭代结果与上一次的迭代结果的最大相差值小于1e-6时认为迭代已经收敛输出结果即可当然也可以换成其它结束迭代方法如判断两个向量之差的二范数。 与直接使用克拉默法则计算准确的解以及matlab计算结果比较不难发现其 x 1 x_1 x1和 x 3 x_3 x3均不为0只是是一个在我们设定的误差范围内接近0的数符合迭代法的解的性质只能在设定的误差范围内得到一个近似的解。
