高斯列主元消去

高斯（Gauss）列主元消去法：对给定的n阶线性方程组A x = b，首先进行列主元消元过程，回代，最后得到方程的解。

代码

C/C++代码：

#include <cmath>
#include <iostream>
using namespace std;
void gauss(double (*a)[5], double *b, double *end, int n);//高斯列主元消去
void swapLine(double (*a)[5], double *b, int n, int begin);//交换矩阵A的行，将每列最大的数所在的一列放在对角线
int main() {
  double a[5][5] = {{0.8147, 0.0975, 0.1576, 0.1419, 0.6557},
                    {0.9058, 0.2785, 0.9706, 0.4218, 0.0357},
                    {0.1270 * (10 ^ 10), 0.5469, 0.9572, 0.9157, 0.8491},
                    {0.9134, 0.9575, 0.4854 * (10 ^ 8), 0.7922, 0.9340},
                    {0.6324, 0.9649, 0.8003, 0.9595, 0.6787}};
    //对两个大数赋值，上面无法赋值
  a[2][0] = 0.1270 * pow(10, 10);
  a[3][2] = 0.4854 * pow(10, 8);
  double b[] = {0.000000002258000, 0.000000001597700, 1.270000002354900,
                0.024270003904200, 0.000000003360250};
  int n = 5;
  double end[5] = {0};
  gauss(a, b, end, n);
  cout << "结果：" << endl;
  for (int i = 0; i < n; i++) {
    cout << "x" << i << "=" << end[i] << "\t" << endl;
  }
}
void gauss(double (*a)[5], double *b, double *end, int n) {
  double x = 0;//中间值
  double y = 0;//中间值
  // 将矩阵化为上三角阵，遍历系数矩阵的每一行

  for (int k = 0, i = 0, j = 0; k < n; k++) {
    swapLine(a, b, n, k);//交换
    // 输出矩阵A
    for (int i = 0; i < n; i++) {
      for (int j = 0; j < n; j++) {
        cout << a[i][j] << " ";
      }
      cout << endl;
    }
    cout << endl;
    // 对从k+1行开始的行进行初等行变换

    y = a[k][k];
    if(y==0){
      return;
    }
    //消元
    for (i = k + 1; i < n; i++) {
        //如果那一列需要消去的第一个已为0,则表示该行不需要消。
      if (a[i][k] == 0) {
        continue;
      }
      x = -y / a[i][k];
      for (j = 0; j < n; j++) {
        a[i][j] = a[i][j] * x + a[k][j];
      }
      b[i] = b[i] * x + b[k];
    }
  }

  // 将上三角矩阵从下到上算出每个x
  for (int i = n - 1, j = 0; i >= 0; i--) {
    x = 0;
    for (j = i + 1; j < n; j++) {
      x += a[i][j] * end[j];
    }
    end[i] = (b[i] - x) / a[i][i];
  }
}
void swapLine(double (*a)[5], double *b, int n, int begin) {
  double max = a[begin][begin];
  int max_index = begin;
  //搜索要交换的行
  for (int j = begin + 1; j < n; j++) {
    // cout<<a[j][begin]<<endl;
    if (abs(max) < abs(a[j][begin])) {
      max = a[j][begin];
      max_index++;
    }
  }

  // cout<<"max_index="<<max_index<<"max:"<<a[max_index][begin]<<endl;
  //如果不需要交换则退出
  if (max_index != begin) {
    for (int i = begin; i < n; i++) {
      max = a[begin][i];
      a[begin][i] = a[max_index][i];
      a[max_index][i] = max;
    }
    max = b[begin];
    b[begin] = b[max_index];
    b[max_index] = max;
  }
}

Python代码：

import numpy as np
import math
from scipy.linalg import solve
a=[[0.8147,0.0975,0.1576,0.1419,0.6557],
    [0.9058 ,0.2785,0.9706 ,0.4218,0.0357],
    [0.1270*(10^10),0.5469 ,0.9572 ,0.9157 ,0.8491],
    [0.9134 ,0.9575 ,0.4854*(10^8) ,0.7922,0.9340],
    [0.6324 ,0.9649 ,0.8003 ,0.9595 ,0.6787]]
a[2][0]=0.1270*math.pow(10,10)
a[3][2]=0.4854*math.pow(10,8)
b=[0.000000002258000,0.000000001597700,1.270000002354900,0.024270003904200,0.000000003360250]
a1=[[4,3,0],[3,4,-1],[0,-1,4]]
b1=[24,30,-24]
end=[0,0,0,0,0]
print(a)
print(-a[0][0]/a[2][0]*a[2][1]+a[0][1])
def gauss(a,b,end,n):
    for k in range(0,n):
        y=a[k][k]
        for i in range(k+1,n):
            if a[i][k]==0:
                continue
            x=-y/a[i][k]
            for j in range(0,n):
                a[i][j]= a[i][j]*x+a[k][j];
            b[i]=b[i]*x+b[k];
    print(a)
    i=n-1
    while i>=0:
        x=0
        j=i+1
        while j<n:
            x+=a[i][j] * end[j]
            j+=1
        end[i]=(b[i]-x) / a[i][i]
        i=i-1
    return end
def swapLine(a,b,begin,n):
    return n
a = np.array(list(a))

#使用scipy库函数来求解
x1 = solve(a, b)
print("scipy lib:",x1)
# sum1=0
# for i in range(0,5):
#     sum1 =sum1 + x1[i]*a[0][i]
# print(sum1)
# print(b[4]/x1[4])
end=gauss(list(a),b,end,5)
print(end)

结果及验证

C/C++代码结果：

结果：
x0=1e-09    
x1=2e-09    
x2=5e-10    
x3=-1e-09   
x4=2e-09

python计算结果：

scipy lib: [ 1.e-09  2.e-09  5.e-10 -1.e-09  2.e-09]
[1e-09, 2.000000000000071e-09, 5.000000000000008e-10, -1.0000000000000489e-09, 1.9999999999999997e-09]

结果一致。