【c++】列主元素高斯消去法求解方程

 

 原代码【截取至大作业4】

void equation_solver(double dF_o[21][21], double F_o[21], double dx[100],int K) {
    int point = 0;
    double temp = 0.0;
    double dF[21][21];double F[21];

    //赋值
    for (int i = 0;i < K;i++) {
        for (int j = 0;j < K;j++) {
            dF[i][j] = dF_o[i][j];
        }
    }
    for (int i = 0;i < K;i++) {
        F[i] = F_o[i];
    }
    //选主元

    for (int k = 0;k < K-1;k++) {
        point = k;
        for (int l = k;l < K;l++) {
            if (fabs(dF[point][k]) < fabs(dF[l][k])) {
                point = l;
            }
        }
        //换行
        temp = 0;
        temp = F[point];
        F[point] = F[k];
        F[k] = temp;
        for (int l = k;l < K;l++) {
            temp = 0;
            temp = dF[point][l];
            dF[point][l] = dF[k][l];
            dF[k][l] = temp;
        }
        //消去过程
        double mid;
        for (int l = k + 1;l < K;l++) {
            mid = dF[l][k] / dF[k][k];
            F[l] = F[l] - mid * F[k];
            for (int p = 0;p < K;p++) {
                dF[l][p] -= mid * dF[k][p];
            }
        }
    }
   /* std::cout << "上三角： " << '\n';
    for (int i = 0;i < 4;i++) {
        for (int j = 0;j < 4;j++) {
            std::cout << "i" << i << "j" << j << "  " << dF[i][j] << "  ";
        }
        std::cout << "\n";
    }*/

    //求解方程
    int M = K - 1;
    dx[M] =F[M] / dF[M][M];
    for (int i = K-2;i > -1;i--) {
        temp = 0;
        //dx[i] = 0;
        for (int j = i + 1;j < K;j++) {
            temp += dF[i][j] * dx[j] / dF[i][i];
            //dx[i] -= dF[i][j] * dx[j] / dF[i][i];
        }
       dx[i] = F[i] / dF[i][i] - temp;
        //std::cout << "dx" << i << ": " << dx[i] << "\n";
    }

}

以书上p17-18 的例1为例。

point1：上述代码中，矩阵的长度是限定的。可以通过std::vector进行改进；

point2：列主元素高斯消去法求解线性方程组需要系数矩阵非奇异。上述代码并未体现；

 

 

#include <iostream>
#include <vector>
#include <algorithm>
#include <math.h>
#include <assert.h>

using vec = std::vector<std::vector<double>>;
using vecRow = std::vector<double>;

void show(vec& A) {
    for (vecRow i : A) {
        for (double j : i) {
            std::cout << j << " ";
        }
        std::cout << '\n';
    }
    /* for (int i = 0;i != A.size();++i) {
         for (int j = 0;j != A[0].size();++j) {
             std::cout << A[i][j] << "  ";
         }
         std::cout << "\n";
     }*/
}

vecRow equationSolver(vec A,vecRow B) {   //如果这里为vec&A vecRow&B ，则A和B的值会改变。如果没有&，则不变，改变的为A和B的副本
    int row1{ static_cast<int> (A.size()) };
    int col1{ static_cast<int>(A[0].size()) };
    int row2{ static_cast<int>(B.size()) };
    assert((row1 == row2) || (row1 > 0) && "请检查输入");
    int point;
    double temp;
    double m;
    vecRow x(row1);
    if (row1 == 1) {
        x[row1] = B[row1] / A[row1][row1];
        return x;
    }
    //一列一列来
    for (int k = 0;k != col1 - 1;++k) {
        //选主元
        point = k;
        //寻找k列中绝对值最大的数的行point
        for (int r = k;r != row1;++r) {
            if (fabs(A[point][k]) < fabs(A[r][k])) {
                point = r;
            }
        }
        //判断A是否为奇异矩阵
        assert(fabs(A[point][k]) > 1e-6 && "A为奇异矩阵");
        //换行
        temp = B[k];
        B[k] = B[point];
        B[point] = temp;
        for (int i = 0;i != col1;++i) {
            temp = A[point][i];
            A[point][i] = A[k][i];
            A[k][i] = temp;
        }
        //消去
        m = 0.0;
        for (int i = k+1;i != col1;++i) {
            m = A[i][k] / A[k][k];
            B[i] -= m * B[k];
            for (int j = 0;j != col1;++j) {
                A[i][j] -= m * A[k][j];
            }
        }
    }
    //解方程
    double sum{ 0.0 };
    x[row1 - 1] = B[row1 - 1] / A[row1 - 1][row1 - 1];
    for (int i = row1-2;i >= 0;--i) {
        sum = 0.0;
        for (int j = i+1;j != col1;++j) {
            sum += A[i][j] * x[j] / A[i][i];
        }
        x[i] = B[i] / A[i][i] - sum;
    }
    /*for (double i : x) {
        std::cout << i << "  ";
    }*/
    return x;
}

int main()
{
    vec A{ {0.012,0.010,0.167},
           {1.000,0.833,5.910},
           {3200 ,1200 ,4.200} };
    vecRow B{ 0.6781,12.1,981 };
    vecRow x;
    x =  equationSolver(A, B);
    for (double i : x) {
        std::cout << i << "  ";
    }
    std::cout << "A:\n";
    show(A);

}