编程最重要的就是要有很强的目的性,面对绝对真诚的计算机时,如果自己还模棱两口,那又如何让计算机明白我们要做的事情呢?所以在编程之前,规划好绝对是磨刀不误砍柴工。笔者基于已有的编程知识,正好又在学习《数值分析》,数学作为编程的心法,重要性不言而喻,所以随着《数值分析》的学习,争取将教材中提到的算法,都用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
第一次写文章,哈哈,接下来随着学习,再进行补充,祝读者每天开开心心。