数值分析（一）：全选主元高斯消元法的C++实现

编程最重要的就是要有很强的目的性，面对绝对真诚的计算机时，如果自己还模棱两口，那又如何让计算机明白我们要做的事情呢？所以在编程之前，规划好绝对是磨刀不误砍柴工。笔者基于已有的编程知识，正好又在学习《数值分析》，数学作为编程的心法，重要性不言而喻，所以随着《数值分析》的学习，争取将教材中提到的算法，都用C++实现一遍，这里实现线性方程的直接解法：全选主元高斯消元法。（完整代码在最后，GAUS类）
（一）、全选高斯消元法的功能和步骤
明确好目的，高斯消元法就是要解线性方程组，获得精确解（实际上还是存在舍入误差），方程如下图：

理清步骤，编程之前先用数学语言将算法的思路表达清楚，每一门学科都是一门特殊的语言，思路表达清楚之后就方便编程实现了。全选主元的高斯消元法步骤用数学语言描述如下：

（二）、每一步的代码实现
1：对于k从0~n-2做如下运算，这是一个大循环；
全选主元，将绝对值最大的元素交换到主对角线上，代码如下：

        d = 0.0;for ( i =k; i <=n-1 ; i++){for ( j = k; j <=n-1; j++){t = a[i][j];if (t > d){d = t;js[k] = j;is = i;}}}if (d == 0.0)flag = false;else{if (js[k] != k){for (i = 0; i < n; i++){t = a[i][k];a[i][k] = a[i][js[k]];a[i][js[k]] = t;}}if (is != k){for ( j = 0; j < n; j++){t = a[k][j];a[k][j] = a[is][j];a[is][j] = t;}t = b[k]; b[k] = b[is]; b[is] = t;}}

然后进行归一化和消元，归一化是将主元以及解向量所在的行全部除以主元素，将对角线上的元素化成1，消元,是将主元素一下的位置全部化成0，主元所在的行以外元素也要跟着变，这里有一点困扰了一段时间，那就是编程和做数学题的区别，实际上这里消元之后，主对角线下边的元素并不是0，计算机压根就没去算它，把它消成0，为什么呢？因为没有意义，后面的回代解x时，完全就没用到对角线以下的位置元素，所以在形式上，对角线下的元素程序并没有变成0，而解数学题的时候，对角线以下的元素还是必须写成0，这也就是为什么消元的时候，下标i，j全是是从k+1开始的，归一化，消元步骤的代码如下所示：

 if (!flag){cout << "该矩阵奇异，不存在唯一解" << endl;delete[] js;return;}d = a[k][k];for ( i = k; i <n; i++){a[k][i] = a[k][i] / d;}b[k] = b[k] / d;for (i = k + 1; i < n; i++){for (j = k+1; j < n; j++){a[i][j] = a[i][j] - a[k][j] * a[i][k];}b[i] = b[i] - b[k] * a[i][k];}

第1步就是将这两段代码写在k从0到n-2的大循环里。
2、进行方程的回代求解，从矩阵下面至上进行回代，从n-1回代到第1行，最终求得所有未知数，数据放在数组b里。

d = a[n - 1][n - 1];
if (d == 0)
{cout << "该矩阵奇异，没有唯一解" << endl;delete[] js;return;
}
b[n - 1] = b[n - 1] / d;
for (i = n - 1; i >= 0; i--)
{t = 0;for (j = i + 1; j < n; j++)t = t + a[i][j] * b[j];b[i] = b[i] - t;
}

3、第3步是将解向量回复成原来的位置，因为在进行主元选择的时候，打乱了解向量中的顺序，所以再调换回来，代码如下：

    js[n - 1] = n - 1;for ( k = n-1; k >=0; k--){if (js[k] != k){t = b[k];b[k] = b[js[k]];b[js[k]] = t;}}

这样，全选主元的高斯消元法就完成了，再把完整的代码附上，头文件GAUS.h,源文件是GAUS.cpp。这里输入矩阵采用gaus.txt文件，输出文件也自己命名，得到计算结果，保存在txt文件里。

    #include <iostream>#include <fstream>#include <cmath>using namespace std;class GAUS
{
private:int n;double **a, *b;public:GAUS(int nn) :n(nn) {a = new double*[n];for (int i = 0; i < n; i++)a[i] = new double[n];b = new double[n];}void input();void gauss();void output();~GAUS() {int i = 0;for (i = 0; i < n; i++)delete[] a[i];delete[] a;delete[] b;}
};//GAUS.cpp
#include "GAUS.h"
void GAUS::input()
{int i = 0, j = 0;char str1[20];cout << "输入文件名" << endl;cin >> str1;ifstream fin(str1);if (!fin)cout << "没有找到该文件" << endl;for (i = 0; i < n; i++){for (j = 0; j < n; j++)fin >> a[i][j];}for (j = 0; j < n; j++)fin >> b[j];fin.close();
}
void GAUS::gauss()
{int i = 0, j = 0, k = 0;int *js, is = 0;js = new int[n];bool flag = true;double t = 0, d = 0;for ( k = 0; k <= n-2; k++){d = 0.0;for ( i =k; i <=n-1 ; i++){for ( j = k; j <=n-1; j++){t = a[i][j];if (t > d){d = t;js[k] = j;is = i;}}}if (d == 0.0)flag = false;else{if (js[k] != k){for (i = 0; i < n; i++){t = a[i][k];a[i][k] = a[i][js[k]];a[i][js[k]] = t;}}if (is != k){for ( j = 0; j < n; j++){t = a[k][j];a[k][j] = a[is][j];a[is][j] = t;}t = b[k]; b[k] = b[is]; b[is] = t;}}if (!flag){cout << "该矩阵奇异，不存在唯一解" << endl;delete[] js;return;}d = a[k][k];for ( i = k; i <n; i++){a[k][i] = a[k][i] / d;}b[k] = b[k] / d;for (i = k + 1; i < n; i++){for (j = k+1; j < n; j++){a[i][j] = a[i][j] - a[k][j] * a[i][k];}b[i] = b[i] - b[k] * a[i][k];}}d = a[n - 1][n - 1];if (d == 0){cout << "该矩阵奇异，没有唯一解" << endl;delete[] js;return;}b[n - 1] = b[n - 1] / d;for (i = n - 1; i >= 0; i--){t = 0;for (j = i + 1; j < n; j++)t = t + a[i][j] * b[j];b[i] = b[i] - t;}js[n - 1] = n - 1;for ( k = n-1; k >=0; k--){if (js[k] != k){t = b[k];b[k] = b[js[k]];b[js[k]] = t;}}delete[] js;
}void GAUS::output()
{int i;char str2[20];cout << "\n输出文件名";cin >> str2;ofstream fout(str2);if (!fout){cout << "/n不能打开文件" << str2 << endl;exit(1);}fout << endl;cout << endl;for ( i = 0; i < n; i++){fout << b[i] << " ";cout << b[i] << " ";}fout << endl;cout << endl;fout.close();
}

笔者主函数调用是：

int main()
{GAUS c(4);c.input();c.gauss();c.output();system("pause");return 0;
}

最后运行一番，方程如下图：

在工程里面建立txt文件，把系数全部写入里面，命名gaus.txt:
0.2368 0.2471 0.2568 1.2671
0.1968 0.2071 1.2168 0.2271
0.1581 1.1675 0.1768 0.1871
1.1161 0.1254 0.1397 0.1490
1.8471 1.7471 1.6471 1.5471
运行程序，输入文件名gaus.txt,输入结果保存的txt文件名，如result.txt,就可以运行解方程了，结果如下：
1.04058 0.987051 0.93504 0.881282
第一次写文章，哈哈，接下来随着学习，再进行补充，祝读者每天开开心心。

数值分析（一）：全选主元高斯消元法的C++实现相关推荐

Matlab实现 LU分解法解线性方程组（全选主元列选主元）
选主元LU分解实验内容:列选主元LU分解和全选主元LU分解求解线性方程组计算方法: 全选主元消元法 1.1 初始化根据参数A.b,记录下矩阵.右端项的尺寸n: 以得到的尺寸n初始化解向量x: 同 ...
全选主元matlab,全选主元高斯-约当法实矩阵求逆及传统方法求逆
/* 算法没有做过什么改动,加上了内存释放的部分. 算法介绍矩阵求逆在3D程序中很常见,主要应用于求Billboard矩阵.按照定义的计算方法乘法运算,严重影响了性能.在需要大量Billboard矩 ...
矩阵求逆（全选主元高斯-约当消去法）
矩阵求逆(实矩阵) 具体的C实现: #include "stdlib.h"#include "math.h"#include "stdio.h&quo ...
数值作业：Guass全选主元消去法之C语言代码
全选主元的基本思想:当变换到第K步时,从系数矩阵的右下角(n-k+1)阶子阵中选取绝对值最大的元素,然后通过行变换与列变换将它交换到主元素a[k][k]的位置上．行变换是不影响最后求解的结果的,但是列 ...
数值计算（一）之解线性方程组（高斯消去法，列选主元消去法，全选主元消去法，杜立特尔分解，克洛特分解，乔里斯基分解）
解线性方程组即解一个多元一次方程组,例如目录消去法分解法消去法原理没有学过高级的解法也没关系,凭借我们初高中的知识足以解决这个问题这是一个多元一次方程组,拥有n个未知量,也有n方程我们 ...
C#，码海拾贝（24）——求解“复系数方程组”的“全选主元高斯－约当消去法”之C#源代码
using System; namespace Zhou.CSharp.Algorithm { /// <summary> /// 求解线性方程组的类 LEquations ...
全选主元高斯消去法 matlab程序,[VB]全选主元高斯消去法求解线性方程组
功能: 使用全选主元高斯消去法求解线性方程组 ' 参数 n - Integer型变量,线性方程组的阶数 ' dblA - Double型 n x n 二维数组,线性方程组的系数矩 ...
全选主元的Gauss-Jordan消元法
选主元的高斯-约当(Gauss-Jordan)消元法在很多地方都会用到,例如求一个矩阵的逆矩阵.解线性方程组(插一句:LM算法求解的一个步骤),等等.它的速度不是最快的,但是它非常稳定(来自网上的定义 ...
全选主元高斯消去法 matlab程序,编译无误的全选主元高斯消去法
算法来源于<> 徐士良著但,我发现他的算法在传递二维数组时,怎么也不能把数据传入.可能是我的水平有限,没有搞明白为什么.但问题总得解决,所以我做了一点点小小的改动.能通过编译,我用 ...

数值分析（一）：全选主元高斯消元法的C++实现

数值分析（一）：全选主元高斯消元法的C++实现相关推荐

最新文章

热门文章