Implement the principle and practice of solving linear equations (matrix, Gaussian elimination)------C + + programming (Advanced article)

Steps:

Where a is a n*n the coefficients of the square vectors x and b are the unknowns and the constant vectors respectively:

The system may have 0, one, or infinitely multiple solutions, depending on the coefficients matrix A and vector B. There are many methods for solving linear systems, where a classical method-Gaussian elimination (https://zh.wikipedia.org/wiki/Gaussian elimination method) is used. First, we exchange A and B to make a into an upper triangular matrix. The upper triangular matrix is 0 for all elements below the diagonal. such as the following form:

It is easy to achieve this goal. To make a (i,j) 0, we multiply it by a constant so that it equals another element on column J, such as Equals a (K,J). Then, with the first equation minus the K equation, a (i,j) becomes 0, and the values of the other elements of the matrix I line change accordingly.

If such a transformation eventually makes all the elements on the diagonal not 0, the equations have a unique solution, which can be obtained by "back-generation". If there is an element of 0, it means that there are 0 or infinitely more solutions to the equation set.

We are now using C + + to represent the above calculation method. First, define the two specific matrix types to use to simplify the program:

typedef numeric_lib::matrix<double, 2>matrix2;   Matrix Library  : http://www.stroustrup.com/Programming/Matrix/Matrix.h  The entire library is defined in the namespace  numeric_lib Medium typedef numeric_lib::matrix<double, 1> Vector;

　　

Next we describe the Gaussian elimination process as a program:

Vector classic_gaussian_elimination (Matrix2 a,vector b) {    classical_elimination (a);    Return Back_substitution (A, b);}

That is, a copy is created for two inputs A and B (using an assignment function), then a function is called to solve the equations, and the return function is called to evaluate the result and return the result. The key point is that the way we decompose the problem and the symbolic representation all come from the original mathematical description. The next thing to do is to implement Classic_elimination () and Back_substitution (), and the solution agreement comes entirely from the math textbook:

void Classical_elimination (matrix2&a,vector& b) {    const Index n=a.dim1 ();    The 1th column has been traversed to the penultimate column    //diagonal So the element is populated with 0 for    (Index j=0;j<n-1;++j) {        Const double pivot =a (J,J);        if (pivot==0) cerr<< "error: There is a diagonal bit of 0" << ' \ n ';        The elements below the diagonal of line I are populated with 0 for         (Index i=j+1;i<n;++i) {        const double mult =a (i,j)/pivot;        A[i].slice (j) =scale_and_add (A[j].slice (j),-mult,a[i].slice (j)); A[i].slice (j) means from a[i][j] to the end of this line.        B (i)-=mult*b (j);//make corresponding changes to B}}}

Pivot indicates that the current line is on a diagonal element, and it must be non-0. Because it needs to be used as a divisor; if it is 0, we discard the calculation and throw an exception:

Vector back_substitution (const matrix2&a,const vector&b) {const Index n=a.dim1 (); Vector x (n); for (Index i=n-1;i>=0;--i) {double s=b (i)-dot_product (A[i].slice (i+1), X.slice (i+1)), if (Double M=a (i,i )) x (i) =s/m;elsethrow back_subst_failure (i);} return x;}

complete Sample program:

#include <iostream> #include "Matrix.h"//matrix Library: http://www.stroustrup.com/Programming/Matrix/Matrix.h# Include "MatrixIO.h"//matrixio Library: Http://www.stroustrup.com/Programming/Matrix/MatrixIO.h provides very simple I/O functionality for one-dimensional two-dimensional using　　namespace Numeric_lib; The entire library is defined in the namespace Numeric_lib using namespace Std;typedef numeric_lib::matrix<double, 2>matrix2; typedef numeric_lib::matrix<double, 1> vector;void classical_elimination (matrix2& A, Vector& b) {const Index n = a.dim1 ();//The 1th column is traversed to the penultimate column//diagonal so the element is filled with 0for (Index j = 0; j<n-1; ++j) {Const double pivot = A (J, J); if (pi Vot = = 0) cerr<< "error: One of the diagonal bits is 0" << ' \ n ';//The elements below the diagonal of line I are populated with 0 for (Index i = j + 1; i<n; ++i) {const double MUL t = A (i, j)/pivot; A[i].slice (j) = Scale_and_add (A[j].slice (j),-mult, A[i].slice (j)), B (i)-= Mult*b (j);//make corresponding changes to B}}}vector back_ Substitution (const matrix2&a, const VECTOR&AMP;B) {Const Index n = a.dim1 (); Vector x (n); for (Index i = n-1; I >= 0; i.) {double s = B (i)-dot_product(A[i].slice (i + 1), X.slice (i + 1)), if (double m = A (i, I)) x (i) = s/m;elsecerr<< "Error: There is a diagonal bit of 0" << ' \ n ';} return x;} Vector Classic_gaussian_elimination (Matrix2 A, vector b) {classical_elimination (A, b); return Back_substitution (A, b);} int main () {double val2[3][3] = {2,1,-1,-3,-1,2,-2,1,2};d ouble val1[3] = {8,-11,-3}; Matrix2 A (VAL2); Vector b (Val1); Cout<<classic_gaussian_elimination (A, b);}

　　

Improved:

The problem with pivot 0 is avoidable, and we can arrange the marching lines so that the 0 and the smaller values are removed from the diagonal, thus obtaining a more robust scheme. "More robust" means not sensitive to rounding errors. However, as we put 0 below the diagonal, the value of the element also changes. Therefore, we have to reorder them over and over again to remove the smaller values from the diagonal (that is, the classic algorithm can be used directly after the matrix is not re-arranged):

void Elim_with_partial_pivot (matrix2& A, vector& b) {const Index n = a.dim1 (); for (Index j = 0; J < N; ++j) {in Dex Pivot_row = j;//Find a suitable principal: for (Index k = j + 1; k < n; ++k) if (ABS (A (k, j)) > abs (A (Pivot_row, j)) Pivot_row = k;//if we find a better principal, swap two lines: if (pivot_row! = j) {a.swap_rows (J, Pivot_row); Std::swap (b (j), B (Pivot_row));}  Elimination: for (Index i = j + 1; i < n; ++i) {const double pivot = A (J, J); if (pivot = = 0) error ("Can ' t solve:pivot==0"); const Double mult = A (i, j)/pivot; A[i].slice (j) = Scale_and_add (A[j].slice (j),-mult, A[i].slice (j)); B (i)-= Mult*b (j);}}

Here we use Swap_rows () and scale_and_multipy (), so that the program is more custom-made, and we don't have to write loop code explicitly.

Random number test:

void Solve_random_system (Index n) {Matrix2 A = Random_matrix (n); Vector B = Random_vector (n), cout << "a=" << A << ' \ n '; cout << "b=" << b << ' \ n '; try {V Ector x = Classic_gaussian_elimination (A, b); cout << "Classic Elim solution is x =" << x << ' \ n '; Vector v = a*x;cout << "a*x=" << v << ' \ n ';} catch (const exception& e) {cerr << e.what () << ' \ n ';}}

The program enters the catch statement in three cases:

• There is a bug in the code.
• The input makes classic_elimination error (Elim_with_partial_pivot can do better in many cases).
• Rounding errors cause problems.

In order to test our program, we output a*x, whose value should be the B-phase order. But given the rounding error, the result is considered correct if the value is close enough to B, which is why the following statement is not used in the test program to determine whether the result is correct:

if (a*b!=b) error ("Substitution failed");

In computers, floating-point numbers are just approximations of real numbers, so we must accept approximate calculations. In general, you should avoid using = = and! = To determine whether it is correct.

Matrix and vector multiplication are not defined in the matrix library, so we define this operation for the test program:

Vector operator* (const matrix2&m, const VECTOR&U) {Const Index n = m.dim1 (); Vector v (n); for (Index i = 0; i < n; ++i) v (i) = Dot_product (M[i], u); return v;}

Random_matrix () and Random_vector () are simple applications of random numbers. Index is an indexed type, which is defined by a typedef.

Complete program:

#include <iostream> #include <random> #include <time.h> #include "Matrix.h"//matrix Library:/http Www.stroustrup.com/Programming/Matrix/Matrix.h#include "MatrixIO.h"//matrixio Library: http://www.stroustrup.com/ Programming/matrix/matrixio.husing namespace numeric_lib;//The entire library is defined in the namespace Numeric_lib using namespace Std;typedef Numeric_lib::matrix<double, 2>matrix2; typedef numeric_lib::matrix<double, 1> vector;void classical_elimination (matrix2& A, Vector& b) {const Index n = a.dim1 ();//The 1th column is traversed to the penultimate column//diagonal so the element is filled with 0for (Index j = 0; j<n-1; ++j) {Const double pivot = A (J, J); if (pi Vot = = 0) cerr<< "error: One of the diagonal bits is 0" << ' \ n ';//The elements below the diagonal of line I are populated with 0 for (Index i = j + 1; i<n; ++i) {const double MUL t = A (i, j)/pivot; A[i].slice (j) = Scale_and_add (A[j].slice (j),-mult, A[i].slice (j)), B (i)-= Mult*b (j);//make corresponding changes to B}}}vector back_ Substitution (const matrix2&a, const VECTOR&AMP;B) {Const Index n = a.dim1 (); Vector x (n); for (Index i = n-1; I >= 0; i.) {double S= B (i)-Dot_product (A[i].slice (i + 1), X.slice (i + 1)), if (double m = A (i, I)) x (i) = S/m;else;} return x;} void Elim_with_partial_pivot (matrix2& A, vector& b) {const Index n = a.dim1 (); for (Index j = 0; J < N; ++j) {in Dex Pivot_row = j;//Find a suitable principal: for (Index k = j + 1; k < n; ++k) if (ABS (A (k, j)) > abs (A (Pivot_row, j)) Pivot_row = k;//if we find a better principal, swap two lines: if (pivot_row! = j) {a.swap_rows (J, Pivot_row); Std::swap (b (j), B (Pivot_row));}  Elimination: for (Index i = j + 1; i < n; ++i) {const double pivot = A (J, J); if (pivot = = 0) error ("Can ' t solve:pivot==0"); const Double mult = A (i, j)/pivot; A[i].slice (j) = Scale_and_add (A[j].slice (j),-mult, A[i].slice (j)); B (i)-= Mult*b (j);}} Vector Classic_gaussian_elimination (Matrix2 A, vector b) {Elim_with_partial_pivot (A, b);//classical_elimination (A, B) ; return Back_substitution (A, b);} Vector operator* (const matrix2&m, const VECTOR&AMP;U) {Const Index n = m.dim1 (); Vector v (n); for (Index i = 0; i < n; ++i) v (i) = Dot_product (M[i], u); return v;} int max0 = 10; Vector Random_vector (Index n) {vector v (n);d efault_random_engine ran{(unsigned int) (time (0) +2)};uniform_int_ Distribution<> ureal{0,max0};for (Index i = 0; i < n; ++i) {V (i) = Ureal (ran);} return v;} Matrix2 Random_matrix (Index N) {Matrix2 V (n,n);d efault_random_engine ran{(unsigned int) time (0)};uniform_int_ Distribution<> ureal{0,max0};for (Index i = 0; i < n; ++i) {for (index j = 0; J < N; ++j) v[i][j] = Ureal (ran );} return v;} void Solve_random_system (Index n) {Matrix2 A = Random_matrix (n); Vector B = Random_vector (n), cout << "a=" << A << ' \ n '; cout << "b=" << b << ' \ n '; try {V Ector x = Classic_gaussian_elimination (A, b); cout << "Classic Elim solution is x =" << x << ' \ n '; Vector v = a*x;cout << "a*x=" << v << ' \ n ';} catch (const exception& e) {cerr << e.what () << ' \ n ';}} int main () {/*double val2[3][3] = {2,1,-1,-3,-1,2,-2,1,2};d ouble val1[3] = {8,-11,-3};Matrix2 A (VAL2); Vector b (Val1); Cout<<classic_gaussian_elimination (A, b); */solve_random_system (4);}

　　

Principles and Practice of C + + programming (Advanced article)

Implement the principle and practice of solving linear equations (matrix, Gaussian elimination)------C + + programming (Advanced article)