©2020-2024 ioDraw All rights reserved
在上一篇博客里面,笔者介绍了解线性方程组的LU分解法,这篇来介绍一个新的方法,迭代法.解线性方程组的迭代法有多种,其中就有Jacobi迭代法,它的原理是什么呢?有如下的线性方程组Ax=b,可将其变形为=>Mx=Nx+b=>x=M
-1Nx+M-1b,设B=M-1N=M-1(M-A)=E-M-1A,f=M-1b,即可得到迭代式:X(k+1)=Bx(k)
+f,这里我们只需要设置一个初始的x向量,依次将前一步的xk代入到迭代式中,就可以得到x(k+1)的结果
关于迭代法的两个注意事项:
1.迭代法相比于其他方法在计算大型稀疏矩阵矩阵方面,是有优势的,但不意味着只能解大型稀疏矩阵
2.并非所有的线性方程组都可以用迭代法进行求解,这是因为不是所有的迭代方程都是收敛的,可能会出现的情况就是,在迭代的过程中,会出现迭代解偏离精确解的情况,并且随着迭代的次数增多,会越偏越大
3.遇到不收敛的情况,就不能用迭代法求解,可以选用前面的Guass消元法或者LU分解法
接下来看代码实现~
老规矩,初始化
double** init_Matrix(int r, int c) { double** p = new double* [r]; int d = c +
1; for (int i = 0; i < r; i++) { p[i] = new double[d]; memset(p[i], 0, sizeof(
double) * d); } cout << "请输入线性方程组对应的增广矩阵:" << endl; for (int i = 0; i < r; i++)
{ for (int j = 0; j < d; j++) { cin >> p[i][j]; } } return p; }
检测是否达到精度要求
bool isRight(double**p,int r,double*x) { double sum1 = 0, flag = 0,sum2=0; for
(int i = 0; i < r; i++) { sum1 = 0,flag=0; for (int j = 0; j < r; j++) { sum1 +=
x[j] * p[i][j]; } flag= fabs(p[i][r] - sum1); if (flag>one_Precision)
//解代入单个方程式的误差过大 { return false; } else { sum2 += flag; } } if (sum2>
total_Precision)//整体误差过大 { return false; } return true; }
这一步就是检测迭代得到的解是否达到了我们所要求的精度,也就是终止条件,这里我设置的检测的标准有两个:一个是将解代入单个方程,计算偏差值,若是大于设定的偏差值,就继续迭代,否则就将其偏差值相加,再次进行判断,若未达到设定的偏差值,就说明已经得到了满足精度要求的解,否则,继续判断,函数里面的精度可以自行调整
开始进行迭代
void Iteration(double**p,int r,double*x,double*xx) { int k = 0;//最大迭代次数 double
sum= 0; while (true) { for (int i = 0; i < r; i++) { sum = 0; for (int j=0;j<r;j
++) { if (j==i) { continue; } else { sum -= p[i][j]* xx[j]; } } x[i] = (p[i][r]
+ sum) / p[i][i]; } for (int i = 0; i < r; i++) { xx[i] = x[i]; } printf(
"第%d次迭代结果为:",++k); for (int i = 0; i < r; i++) { printf("%f\t", x[i]); } cout <<
endl; if (k>=MAX_time) { cout << "超出迭代次数上限!停止迭代" << endl; return; } if (isRight
(p, r, x))//精度符合要求 { cout << "精度符合要求,停止迭代,共迭代:" << k << "次" << endl; return; } }
}
每次迭代打印其解向量的值,然后进行精度判断,若不符合要求,则继续迭代,同时为了防止因出现不收敛的情况而导致死循环的情况,设置了一个次数为300的最大迭代次数
综合运用
void Jacobi_main() { int i = 0, j = 0; cout << "请输入线性方程组对应系数矩阵的行和列:" << endl;
cin>> i >> j; double** p = init_Matrix(i, j); double* X = new double[i];
//第n+1次跌代 double* x = new double[i];//第n次迭代 memset(x, 0, sizeof(double) * i);
memset(X, 0, sizeof(double) * i); Iteration(p, i, X,x); for (int i = 0; i < j; i
++) { delete[]p[i]; } delete []p; delete []x; }
这里动态分配两个数组,一个存储第n+1次迭代的解,一个存储第n次迭代的解,同时在计算结束后,记得释放动态分配的内存,防止内存泄漏
完整代码及测试数据
#include<iostream> #include<cmath> #include<Windows.h> #define one_Precision
1e-5//单个方程误差 #define total_Precision 3e-5//整体方程误差 #define MAX_time 1000//最大迭代次数
using namespace std; /* 测试数据 3 3 10 -1 0 9 -1 10 -2 7 0 -2 10 6 3 3 20 -1 2 74
2 8 1 -4 1 -2 4 56 3 3 8 -3 2 20 4 11 -1 33 2 1 4 12 5 5 28 -3 0 0 0 10 -3 38
-10 0 -5 0 0 -10 25 -15 0 0 0 0 -15 45 0 0 0 -5 0 0 30 0 */ double** init_Matrix
(int r, int c) { double** p = new double* [r]; int d = c + 1; for (int i = 0; i
< r; i++) { p[i] = new double[d]; memset(p[i], 0, sizeof(double) * d); } cout <<
"请输入线性方程组对应的增广矩阵:" << endl; for (int i = 0; i < r; i++) { for (int j = 0; j < d;
j++) { cin >> p[i][j]; } } return p; } //检测是否为精确解 设置合格精度为10的-3方 bool isRight(
double**p,int r,double*x) { double sum1 = 0, flag = 0,sum2=0; for (int i = 0; i
< r; i++) { sum1 = 0,flag=0; for (int j = 0; j < r; j++) { sum1 += x[j] * p[i][j
]; } flag= fabs(p[i][r] - sum1); if (flag> one_Precision)//解代入单个方程式的误差过大 {
return false; } else { sum2 += flag; } } if (sum2> total_Precision)//整体误差过大 {
return false; } return true; } //用一维数组x存储每次迭代的解向量 void Iteration(double**p,int r
,double*x,double*xx) { int k = 0;//最大迭代次数 double sum = 0; while (true) { for (
int i = 0; i < r; i++) { sum = 0; for (int j=0;j<r;j++) { if (j==i) { continue;
} else { sum -= p[i][j]* xx[j]; } } x[i] = (p[i][r] + sum) / p[i][i]; } for (int
i= 0; i < r; i++) { xx[i] = x[i]; } printf("第%d次迭代结果为:",++k); for (int i = 0; i
< r; i++) { printf("%f\t", x[i]); } cout << endl; if (k>= MAX_time) { cout <<
"超出迭代次数上限!停止迭代" << endl; return; } if (isRight(p, r, x))//精度符合要求 { cout <<
"精度符合要求,停止迭代,共迭代:" << k << "次" << endl; return; } } } void Jacobi_main() { int i
= 0, j = 0; cout << "请输入线性方程组对应系数矩阵的行和列:" << endl; cin >> i >> j; double** p =
init_Matrix(i, j); double* X = new double[i];//第n+1次跌代 double* x = new double[i]
;//第n次迭代 memset(x, 0, sizeof(double) * i); memset(X, 0, sizeof(double) * i);
Iteration(p, i, X,x); for (int i = 0; i < j; i++) { delete[]p[i]; } delete []p;
delete []x; } int main(void) { Jacobi_main(); system("pause"); return 0; }