解线性方程组迭代法之Jacobi迭代法及其算法实现
在上一篇博客里面,笔者介绍了解线性方程组的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 > 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 (1e-5))//解代入单个方程式的误差过大
{
return false;
}
else
{
sum2 += flag;
}
}
if (sum2>(1e-4))//整体误差过大
{
return false;
}
return true;
}
这一步就是检测迭代得到的解是否达到了我们所要求的精度,也就是终止条件,这里我设置的检测的标准有两个:一个是将解代入单个方程,计算偏差值,若是大于设定的偏差值,就继续迭代,否则就将其偏差值相加,再次进行判断,若未达到设定的偏差值,就说明已经得到了满足精度要求的解,否则,继续判断,函数里面的精度可以自行调整
开始进行迭代
void Iteration(double**p,int r,double*x,double*xx)
{
int k = 0,max_time=300;//最大迭代次数
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 <=max_time)
{
cout << "超出迭代次数上限!停止迭代" << endl;
return;
}
if (isRight(p, r, x))//精度符合要求
{
cout << "精度符合要求,停止迭代,共迭代:" << k << "次" << endl;
return;
}
}
}
每次迭代打印其解向量的值,然后进行精度判断,若不符合要求,则继续迭代,同时为了防止因出现不收敛的情况而导致死循环的情况,设置了一个次数为300的最大迭代次数
综合运用
void Jacobi_main()
{
int i = 0, j = 0;
cout << "请输入线性方程组对应系数矩阵的行和列:" <> 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
#include
#include
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 > 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 (1e-5))//解代入单个方程式的误差过大
{
return false;
}
else
{
sum2 += flag;
}
}
if (sum2>(3e-5))//整体误差过大
{
return false;
}
return true;
}
//用一维数组x存储每次迭代的解向量
void Iteration(double**p,int r,double*x,double*xx)
{
int k = 0,max_time=300;//最大迭代次数
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 <=max_time)
{
cout << "超出迭代次数上限!停止迭代" << endl;
return;
}
if (isRight(p, r, x))//精度符合要求
{
cout << "精度符合要求,停止迭代,共迭代:" << k << "次" << endl;
return;
}
}
}
void Jacobi_main()
{
int i = 0, j = 0;
cout << "请输入线性方程组对应系数矩阵的行和列:" <> 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;
}
作者：爱你是长久之计~
