Gauss Jordan elimination method, column principal element elimination method


Gauss elimination method for solving linear equations:

The corresponding elementary row transformation is made for the augmented matrix. After the coefficient matrix on the right is transformed into the unit matrix, the rightmost column is the solution of the original linear equations.

Example: P3389 [template] Gauss elimination method

1. Gauss Jordan (G-J) elimination method

Direct elimination to diagonal matrix
(1) Starting from the first line, find the principal element in sequence (cannot be 0), and move the principal element to the main diagonal through line transformation
(2) Divide all elements in the row where the main element is located by the main element to make the main element into 1 (or after the fourth part)
(3) Other rows subtract the row of the primary element and multiply it by a certain multiple, so that other elements in the column of the primary element other than the primary element are 0
(4) Find the next principal element in the remaining matrix range and repeat the above steps

#include<algorithm>
#include<iostream>
#include<cstdio>
#include<cstring>
#include<cmath>
#include<vector>
using namespace std;
//#define int long long
#define ll long long
#define uql unsigned long long
#define pii pair<int,int>
#define mid ((l + r)>>1)
#define chl (root<<1)
#define chr (root<<1|1)
#define lowbit(x) ( x&(-x) )
const int manx = 100+10;
const int manx2 = 4e7 + 10;
const ll INF = 1e18;
const ll mod = 1000000007;
const double eps=1e-10;
//The loss of precision during the operation results in that the machine representation of the value that should have been 0 is not 0, so ABS (x) < EPS is used to judge 0
int n;
double f[manx][manx];
int Gauss()
{
    for(int i=1;i<=n;i++){
        for(int j=i;j<=n;j++){/********Corresponding step (1)********/
            if(fabs(f[j][i])>eps){
                swap(f[i],f[j]);
                break;
            }
        }
        if(fabs(f[i][i])<eps){//If the columns are all 0, there is no solution
            puts("No Solution");
            return 0;
        }
        for(int j=1;j<=n;j++){/********Corresponding step (3)********/
            if(i==j)continue;
            double temp=f[j][i]/f[i][i];
            for(int k=i;k<=n+1;k++)
                f[j][k]-=temp*f[i][k];
        }
    }
    for(int i=1;i<=n;i++)//Step (2) can also finally convert the number on the main diagonal to 1
        f[i][n+1]/=f[i][i];
    return 1;
}
int main()
{
    scanf("%d",&n);
    for(int i=1;i<=n;i++)
        for(int j=1;j<=n+1;j++)
            scanf("%lf",&f[i][j]);
    if(!Gauss())return 0;
    for(int i=1;i<=n;i++)
        printf("%.2f\n",f[i][n+1]);
    return 0;
}

2. Column principal component elimination method

First eliminate the upper triangular matrix, and then bring it back for solution
(1) Starting from the first column row, select the element with the largest absolute value in the current column as the main element of this column, and move the main element to the main diagonal through row transformation
(2) Divide all elements in the row where the main element is located by the main element to make the main element into 1 (or after the fourth part)
(3) Subtract the row of the principal element from the other rows in the remaining matrix range and multiply it by a certain multiple, so that the elements other than the principal element in the column of the remaining matrix principal element are 0
(4) Find the next principal element in the remaining matrix range and repeat the above steps
(5) Back band solution

#include<algorithm>
#include<iostream>
#include<cstdio>
#include<cstring>
#include<cmath>
#include<vector>
using namespace std;
//#define int long long
#define ll long long
#define uql unsigned long long
#define pii pair<int,int>
#define mid ((l + r)>>1)
#define chl (root<<1)
#define chr (root<<1|1)
#define lowbit(x) ( x&(-x) )
const int manx = 100+10;
const int manx2 = 4e7 + 10;
const ll INF = 1e18;
const ll mod = 1000000007;
const double eps=1e-18;
int n;
double f[manx][manx];
int Gauss()
{
    int Max;
    for(int i=1;i<=n;i++){
        Max=i;
        for(int j=i+1;j<=n;j++){/********Corresponding step (1)********/
            if(f[Max][i]>fabs(f[j][i]))
                Max=j;
        }
        if(fabs(f[Max][i])<eps){
            puts("No Solution");
            return 0;
        }
        swap(f[Max],f[i]);
        for(int j=i+1;j<=n;j++){/********Corresponding step (3)********/
            double temp=f[j][i]/f[i][i];
            for(int k=i;k<=n+1;k++)
                f[j][k]-=temp*f[i][k];
        }
    }
    for(int i=n;i>=1;i--){
        for(int j=i+1;j<=n;j++)/********The corresponding step (5) brings back the solution********/
            f[i][n+1]-=f[i][j]*f[j][n+1];
        f[i][n+1]/=f[i][i];/********Corresponding step (2)********/
    }
    return 1;
}
int main()
{
    scanf("%d",&n);
    for(int i=1;i<=n;i++)
        for(int j=1;j<=n+1;j++)
            scanf("%lf",&f[i][j]);
    if(!Gauss())return 0;
    for(int i=1;i<=n;i++)
        printf("%.2f\n",f[i][n+1]);
    return 0;
}

3. Comparison of several elimination methods

  • Gauss Jordan elimination method: there is no loopback process and the code is simple; Easy to understand. However, because other rows have to be eliminated every time, the efficiency is lower than that of column principal elements; In addition, the absolute value of the main diagonal element will be smaller than that of the elements below the diagonal in the same column when looking for the principal element in sequence, resulting in the extreme phenomenon that the number with the smallest absolute value is used as the denominator in the elimination process, so that the error in the solution process continues to accumulate and expand (therefore, a larger number can also be selected here as the principal element to reduce the error).
  • Column principal component elimination method: one more carry back process. But the efficiency and accuracy are higher than Gauss Jordan elimination method
  • All selected principal component elimination method : the accuracy is higher than that of the column principal element elimination method, but the column transformation is also carried out for the matrix by selecting all principal elements. The column transformation makes the position transformation of the unknowns during the back substitution. It is necessary to transform the order of the Unknowns at the same time while performing the column transformation, which increases the amount of calculation, and the actual time is close to twice that of the column principal element.

4. Accuracy of double

(1) Storage method of double:

(2) Valid bits: 16 ~ 15 bits
Because the last digit has rounding error, it can ensure that at least the first 15 digits are recorded correctly:
(according to the storage format, the base number of 64 bit floating-point number has 52 bits, plus one hidden bit, a total of 53 bits. The maximum valid value of 64 bit floating-point number is 253 = 9007199254740992)

Decimal number to binary number: conversion of integer part: divide by 2 and take remainder; Conversion of decimal part: multiply by 2 for remainder

Sometimes, in the conversion, because the binary digits of the integer part have been greater than 53, or some bits of the binary decimal will repeat repeatedly, resulting in infinity. However, the representation of the computer is limited, so only a certain accuracy can be intercepted in the computer
(error size: The significant number of a floating-point number)

(3) Influence of absolute value of principal component on accuracy in Gaussian elimination

If a smaller number is selected as the principal element, the coefficient temp multiplied by other lines may be too large. At this time, the f[i][k] stored in the computer has error. The greater the temp multiplied by it, because the greater the error accumulated by f[i][k]

reference resources:
1,How to solve linear equations by Gauss elimination method
2,Gauss column principal component elimination method for solving linear equations
3,Column principal component elimination method and total principal component elimination method
4,Gaussian elimination with different principal components
5,What is the accuracy of Gaussian elimination?
6,Storage mode of floating point number in computer

Added by Mactek on Tue, 08 Mar 2022 14:19:05 +0200