数值分析(一):全选主元高斯消元法的C++实现

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

(二)、每一步的代码实现
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;
			}
		}

免费领qq红包活动然后进行归一化和消元,归一化是将主元以及解向量所在的行全部除以主元素,将对角线上的元素化成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 

第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 
    #include 
    #include 
    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 = 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;
}

最后运行一番,方程如下图:
图3 待解线性方程组
在工程里面建立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
第一次写文章,哈哈,接下来随着学习,再进行补充,祝读者每天开开心心。

猜你喜欢

转载自blog.csdn.net/baidu_41647951/article/details/83010889