In the previous blog, the author introduced LU decomposition method for solving linear equations. This paper introduces a new method, iterative method. There are many iterative methods for solving linear equations, including Jacobi iterative method. What is its principle? There are the following linear equations, Ax=b, which can be transformed into = > MX = NX + B = > x = m-1nx + m-1b. Let B=M-1N=M-1(M-A)=E-M-1A,f=M-1b, Then we can get the iterative formula: X(k+1)=Bx(k)+f. here we just need to set an initial x vector, and then substitute xk of the previous step into the iterative formula in turn, then we can get the result of x(k+1)
There are two considerations about the iterative method:
1. Compared with other methods, iterative method has advantages in calculating large sparse matrix, but it does not mean that only large sparse matrix can be solved
2. Not all linear equations can be solved by iterative method, because not all iterative equations are convergent, and the possible situation is that in the process of iteration, the iterative solution deviates from the exact solution, and as the number of iterations increases, the deviation will increase
3. In case of non convergence, it can't be solved by iterative method. You can choose the former Guass elimination method or LU decomposition method
Let's see the code implementation~
Old rules, initialization
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 << "Please input the corresponding augmented matrix of linear equations:" << endl; for (int i = 0; i < r; i++) { for (int j = 0; j < d; j++) { cin >> p[i][j]; } } return p; }
Check whether the accuracy requirements are met
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>(1e-5))//Too much error in solving single equation { return false; } else { sum2 += flag; } } if (sum2>(1e-4))//Overall error is too large { return false; } return true; }
This step is to check whether the solution obtained by iteration meets the precision required by us, that is, the termination condition. Here, I set two standards for detection: one is to substitute the solution into a single equation and calculate the deviation value. If it is greater than the set deviation value, continue the iteration, otherwise add the deviation value and judge again. If it fails to meet the set deviation value, It means that the solution meeting the accuracy requirements has been obtained. Otherwise, continue to judge, and the accuracy in the function can be adjusted by itself
Start iteration
void Iteration(double**p,int r,double*x,double*xx) { int k = 0,max_time=300;//Maximum number of iterations 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("The first%d The result of the second iteration is:",++k); for (int i = 0; i < r; i++) { printf("%f\t", x[i]); } cout << endl; if (k>=max_time) { cout << "The maximum number of iterations is exceeded!Stop iteration" << endl; return; } if (isRight(p, r, x))//The accuracy meets the requirements { cout << "The accuracy meets the requirements,Stop iteration,Co iteration:" << k << "second" << endl; return; } } }
Each iteration prints the value of its solution vector, and then makes the accuracy judgment. If it does not meet the requirements, the iteration will continue. At the same time, in order to prevent the situation of dead cycle caused by non convergence, a maximum number of iterations of 300 is set
Comprehensive application
void Jacobi_main() { int i = 0, j = 0; cout << "Please enter the row and column of the coefficient matrix corresponding to the system of linear equations:" << endl; cin >> i >> j; double** p = init_Matrix(i, j); double* X = new double[i];//The n+1 decline double* x = new double[i];//Nth iteration 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; }
Two arrays are dynamically allocated here, one stores the solution of the nth + 1st iteration and the other stores the solution of the nth iteration. At the same time, remember to release the dynamically allocated memory after the calculation to prevent memory leakage
Complete code and test data
#include<iostream> #include<cmath> #include<Windows.h> using namespace std; /* test data 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 << "Please input the corresponding augmented matrix of linear equations:" << endl; for (int i = 0; i < r; i++) { for (int j = 0; j < d; j++) { cin >> p[i][j]; } } return p; } //Check whether to set - 3 square with qualified precision of 10 for accurate solution 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>(1e-5))//Too much error in solving single equation { return false; } else { sum2 += flag; } } if (sum2>(3e-5))//Overall error is too large { return false; } return true; } //Using one dimensional array x to store the solution vector of each iteration void Iteration(double**p,int r,double*x,double*xx) { int k = 0,max_time=300;//Maximum number of iterations 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("The first%d The result of the second iteration is:",++k); for (int i = 0; i < r; i++) { printf("%f\t", x[i]); } cout << endl; if (k>=max_time) { cout << "The maximum number of iterations is exceeded!Stop iteration" << endl; return; } if (isRight(p, r, x))//The accuracy meets the requirements { cout << "The accuracy meets the requirements,Stop iteration,Co iteration:" << k << "second" << endl; return; } } } void Jacobi_main() { int i = 0, j = 0; cout << "Please enter the row and column of the coefficient matrix corresponding to the system of linear equations:" << endl; cin >> i >> j; double** p = init_Matrix(i, j); double* X = new double[i];//The n+1 decline double* x = new double[i];//Nth iteration 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; }
