c++面向对象的方式打开线性方程组求解：Jacobi与Gausss-Seidel迭代法、高斯消元法

学校课程要求

当然上述代码各大网站已经漫天飞了，随便搜几个回来自己整合一下轻松就完成了，放几个我参考的文章：

雅各比迭代法： https://blog.csdn.net/xiaowei_cqu/article/details/8585703

Gausss-Seidel迭代法： https://blog.csdn.net/qq_27508477/article/details/83152993

高斯消元法： https://blog.csdn.net/qq_26025363/article/details/53027926

但是整合起来的代码风格各异，c风格的c++代码、各种幻数、冗余的表达等等全部挤在一个cpp文件里，看得人眼花缭乱、头晕目眩。于是借此机会正好可以熟悉练习一下之前阅读《c++程序设计语言》----Bjarne Stroustrup,4th edition 学到的面向对象的程序排布和模块化编程，正好最近也在做一个项目需要用到这些好的方法。

我先定义了两个基类 IterativeMethod 和 GaussEli 分别用来存放两种迭代法（Jacobi和Gauss-Seidel）和高斯消元法的算法代码并留出接口用来输入矩阵和精度参数，然后在main函数的文件定义了类 deploy 来进行方程的参数输入、存放和方程解的调用，类 deploy继承自上述两个类 IterativeMethod 和 GaussEli ，所以可以直接调用基类的算法，从自己变量中储存的数据中输入算法所需参数，并输出结果。

下面上c++代码：

//"FieldPractice.cpp"
//This file includes function main, the program starts here
#include"GaussEli.h"
#include"IterativeMethod.h"
#include<cstdlib>
#include<vector>
//Class that deploys all the methods
class deploy:IterativeMethod,GaussEli {
public:
    //Set up basic matrix
    deploy(double deli):delicacy(deli) {
        cout << "input the number of unknows n:";
        cin >> n;
        cout << endl;
        vector<vector<double> > A1(n, vector<double>(n, 0));
        A = A1;
        vector<double> B1(n, 0);
        B = B1;
        cout << "input the matrix:" << endl;
        for (int i = 0; i < n; i++) {
            for (int j = 0; j < n; j++) {
                cin >>A[i][j];
            }
        }
        cout << endl;
        cout << "input b vector:" << endl;
        for (int k = 0; k < n; k++) {
            cin >> B[k];
        }
        cout << endl;
        cout << "formula to break:" << endl;
        for (int a = 0; a < n; a++) {
            for (int b = 0; b < n; b++) {
                cout << A[a][b] << " ";
            }
            cout << "    " << B[a] << endl;
        }
        cout << endl;
    };

    void jacobi() {
        Jacobi(A, B, n, delicacy);
    }
    void gauss() {
        Gauss(A, B, n, delicacy);
    }
    void gauss_direct() {
        Gaussin(A, B, n, delicacy);
    }
private:
    int n;
    double delicacy;
    vector<vector<double> > A;
    vector<double> B;
};

int main() {
    //Input accuracy factor
    deploy dm(pow(10,-7));
    //Iterative method
    dm.jacobi();
    dm.gauss();
    //Direct method
    dm.gauss_direct();
    return 0;
}

//"GaussEli.h"
//This file contains the declaration of the class DirectMethod
#pragma once
#include<vector>
#include<iostream>
#include<cmath>
#define Length 10
using namespace std;
class GaussEli {
public:
    GaussEli() {};
    ~GaussEli() {};
    void Gaussin(vector<vector<double> > A, vector<double> B, int n, double epslong);
private:

};

#pragma once
//"GaussEli.cpp"
//This file contains the definition of "GaussEli.h"
#include"GaussEli.h"
#include<iostream>
using namespace std;

void GaussEli::Gaussin(vector<vector<double> > a, vector<double> b, int n, double epslong)
{
    cout << endl << endl << "Gaussin Elimination:" << endl;
    //判断能否用高斯消元法，如果矩阵主对角线上有0元素存在是不能用的
    for (int i = 0; i < n; i++)
        if (a[i][i] == 0)
        {
            cout << "can't use gaussin method" << endl;
            return;
        }

    int i, j, k;
    vector<double> c(n,0);    //存储初等行变换的系数，用于行的相减
    //消元的整个过程如下，总共n-1次消元过程。
    for (k = 0; k < n - 1; k++)
    {
        //求出第K次初等行变换的系数
        for (i = k + 1; i < n; i++)
            c[i] = a[i][k] / a[k][k];

        //第K次的消元计算
        for (i = k + 1; i < n; i++)
        {
            for (j = 0; j < n; j++)
            {
                a[i][j] = a[i][j] - c[i] * a[k][j];
            }
            b[i] = b[i] - c[i] * b[k];
        }
    }

    //解的存储数组
    vector<double> x(n, 0);
    //先计算出最后一个未知数；
    x[n - 1] = b[n - 1] / a[n - 1][n - 1];
    //求出每个未知数的值
    for (i = n - 2; i >= 0; i--)
    {
        double sum = 0;
        for (j = i + 1; j < n; j++)
        {
            sum += a[i][j] * x[j];
        }
        x[i] = (b[i] - sum) / a[i][i];
    }

    cout << " the solution of the equations is:" << endl;
    cout << endl;
    for (i = 0; i < n; i++)
        cout << "x" << i + 1 << "=" << x[i] << endl;

}

#pragma once
//"IterativeMethod.h"
//This file contains the declaration of the class IterativeMethod
#include<iomanip>
#include<string>
#include<vector>
#include<iostream>
#include<cmath>
using namespace std;
class IterativeMethod {
public:
    IterativeMethod() {};
    ~IterativeMethod() {};
    void Jacobi(vector<vector<double> > A, vector<double> B, int n, double epslong);
    void Gauss(vector<vector<double> > A, vector<double> B, int n, double epslong);
private:
    double MaxOfList(vector<double> x);
};

//"IterativeMethod.cpp"
//This file contains the definition of "ItearativeMethod.h"
#include"IterativeMethod.h"

//函数求数组中的最大值
double IterativeMethod::MaxOfList(vector<double>x) {
    double max = x[0];
    int n = x.size();
    for (int i = 0; i < n; i++)
        if (x[i] > max) max = x[i];
    return max;
}

//雅可比迭代公式
void IterativeMethod::Jacobi(vector<vector<double> > A, vector<double> B, int n, double epslong) {
    vector<double> X(n, 0);
    vector<double> Y(n, 0);
    vector<double> D(n, 0);
    cout << "Jacobi:" << endl;
    int k = 0; //记录循环次数
    do {
        X = Y;
        for (int i = 0; i < n; i++) {
            double tem = 0;
            for (int j = 0; j < n; j++) {
                if (i != j) tem += A[i][j] * X[j];
            }
            Y[i] = (B[i] - tem) / A[i][i];
            cout << left << setw(8) << Y[i] << " ";
        }
        cout << endl;
        k++;
        if (k > 100) {
            cout << "Failure" << endl;
            return;
        }

        for (int a = 0; a < n; a++) {
            D[a] = X[a] - Y[a];
        }
    } while (MaxOfList(D) > epslong || MaxOfList(D) < -epslong);

    return;
}
//Gauss-Sediel method
void IterativeMethod::Gauss(vector<vector<double> > A, vector<double> B, int n, double epslong) {
    vector<double> X(n, 0);
    vector<double> X1(n, 0);
    vector<double> Y(n, 0);
    vector<double> D(n, 0);
    cout << "Gauss:" << endl;
    int k = 0; //记录循环次数
    do {
        X1 = Y;
        for (int i = 0; i < n; i++) {
            double tem = 0;
            for (int j = 0; j < n; j++) {
                if (i != j) tem += A[i][j] * X[j];
            }
            Y[i] = (B[i] - tem) / A[i][i];
            X = Y;
            cout << left << setw(8) << Y[i] << " ";
        }
        cout << endl;
        k++;
        if (k > 100) {
            cout << "Failure!" << endl;
            return;
        }

        for (int a = 0; a < n; a++) {
            D[a] = X1[a] - Y[a];
        }
    } while (MaxOfList(D) > epslong || MaxOfList(D) < -epslong);

    return;
}

 

 

这样的写法比起原来有几个优点：一是分离了算法和方程矩阵的储存，如果后续想加入新的算法也不用重新写一遍数据输入或者重新输入数据，直接用deploy类再调用新的算法类即可；二是将方程矩阵相关的数据存储和输入全部封装到一个类型中，并且输入的过程直接写在 deploy 类的构造函数，遵循了RAII原则（资源获取即初始化），可以避免一些潜在的错误发生；三是如果程序想要求解多个方程而不是一个，可以直接生成新的 deploy 实体即可，避免了按照代码原来的写法还要繁琐的定义一大堆变量的情况；还有模块化编程在debug时候的便利性我就不在此赘述了。

其实上述程序还有可以进一步改进的空间，应该在类 deploy 中增加变量以储存方程的解，并通过调用算法函数的返回值获得解，再进一步进行输出。换句话说，程序的输出即方程的解可以储存在派生类中，而不是代表算法的基类，算法应该只进行计算而其他部分应该留给别的模块来执行，这样更能体现模块化编程的思想。不过这也只是课程作业的一个随笔，顺手记录，余下的部分就交给有兴趣的读者自行完成了。