本次采用全主元求解线性方程组,比上次的列主元消去法http://blog.csdn.net/qq_26025363/article/details/53044843更加的精确,上次列主元只是选出列中最大的那个数,这次选出行、列中最大的那个数当作主元,无疑增加了可靠性。
本次的算法是完全不同的算法,本次考虑了很多东西,节省了很多空间。输入n*n阶矩阵,和n阶向量,本次算法会改变输入矩阵和输出矩阵的值,意味着消元的过程直接是立地消元,该算法节省了很多空间。唯一的空间消耗是动态分配了n维数组,存放交换的行的值。并且本算法改善了判断的过程,如果系数矩阵奇异,就不能解,就返回值返回-1,如果能解,就返回1。可以根据返回值来判断能否解。其次
,如果还有待改善的话,
1.我认为可以用抛出异常来代替处理不能解的情况
2.可以看到,我用来存放交换的行的动态数组与col的定义是重的,可以删去col,但是对于col这个只占4个字节的内存,影响不大,用来反而能增加程序可读性。
下面贴代码:
//全选主元高斯消去法求解实系数线性方程组
#include<iostream>
#include<cmath>
using namespace std;
const int n=4; //定义矩阵的阶数
//交换元素
void SWAP(double& a, double& b )
{
double c;
c = a;
a = b;
b = c;
}
//由于矩阵运算会有乘除法,所以不能用int数组存储数,选择double
int Rgauss(double a[n][n], double b[n])
{
int i, j, k, l = 11; //i、j、k为循环变量,l是个flag,标志是否是奇异矩阵
int* ex = new int[n]; //用于记忆交换的数组,为了解向量顺序输出 //new一个数组存放
int col, row; //存放第k次选出主元的行和列
double maxs; //存放系数元素的绝对值
for (k = 0; k <= n - 2; k++) //进行k-2次循环消元
{
//选最大元素的过程
double max = 0;
ex[k] = k;
for (i = k; i <= n - 1; i++)
for (j = k; j <= n - 1; j++)
{
maxs = fabs(a[i][j]);
if (maxs > max)
{
max = maxs;
ex[k] = j;
row = i;
col = j;
}
}
//判断是否奇异,非奇异可以继续执行
if (max = 0) l = 0;
else
{
if (col != k) //交换列
for (i = 0; i < n; i++)
SWAP(a[i][k], a[i][col]);
if (row != k) //交换行
for (i = k; i < n; i++)
{
SWAP(a[k][i], a[row][i]);
SWAP(b[k], b[row]);
}
}
//如果奇异就返回-1,代表失败,并释放内存
if (l == 0)
{
delete[] ex;
cout << "failure" << endl;
return -1;
}
//消元过程
double s;
for (j = k + 1; j < n; j++)
{
s = a[j][k] / a[k][k];
cout << "s " << s << endl;
for (i = k+1; i < n; i++)
a[j][i] = a[j][i] - s*a[k][i];
b[j] = b[j] - s*b[k];
}
}
// 如果最后一个元素是0,则方程组无解。
if (0 == a[n - 1][n - 1])
{
delete[] ex;
cout << "failure" << endl;
return -1;
}
//回带求解过程
b[n - 1] = b[n - 1] / a[n - 1][n - 1];
cout << "b[n-1] " << b[n-1] << endl;
for (i = n - 2; i >= 0; i--)
{
double sum = 0.0;
for (j = n - 1; j > i; j--)
{
sum += a[i][j] * b[j];
}
b[i] = (b[i] - sum) / a[i][i];
cout << "sum= " << sum;
cout << endl;
}
//恢复解的正常次序
int xx;
ex[n - 1] = n - 1;
for (k = n - 1; k >= 0; k--)
{
if (ex[k] != k)
{
xx = ex[k];
SWAP(b[k], b[xx]);
}
}
delete[] ex; //释放内存
return 1; //返回成功的值
}
int main()
{
double a[4][4] = { 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 };
double b[4] = { 1.8471,1.7471,1.6471,1.5471 };
//结果是1.04058,0.987051,0.93504,0.881282
//double a[3][3] = { 1,2,2,4,4,2,4,6,4 };
//double b[3] = { 3,6,10 };
//测试结果是-1,3,-1;
if(Rgauss(a, b) != -1) //判断是否能解,能解就输出
{
for (int i = 0; i < n; i++)
cout << b[i] << endl;
}
else
cout << "can't use this way" << endl;
return 0;
}
如果觉得本文有问题,欢迎探讨!!