点击这里可以跳转至
【1】矩阵汇总：http://www.cnblogs.com/HongYi-Liang/p/7287369.html
【2】矩阵生成：http://www.cnblogs.com/HongYi-Liang/p/7275278.html
【3】矩阵加减：http://www.cnblogs.com/HongYi-Liang/p/7287403.html
【4】矩阵点乘：http://www.cnblogs.com/HongYi-Liang/p/7287324.html
【5】矩阵化简：现在的位置
（待续）
...
C++语言：
高斯消元法：
继续使用这个矩阵
当我们使用高斯消元（无回代）化简这个矩阵，是这样算的：
上述过程归纳为：
- 找到第一行行的主元（第一行第一个数：1）
- 消除第而三行的的第一个数（r2-2*r1;r3-4*r1）
- 找到第二行的主元（第二行第二个数：-2）
- 消除第三行的第二个数（r3-3/2*r2）
可以发现实际上是1和2两个步骤的循环，所以写成循环的形式
- 从第一行开始到最后一行
- 找主元：找出第i的主元(第i行第i个数)
- 消元：消除下面所有行的第i个数（下面每一行减去x倍的第一行来消除第i列）
到目前为止，基本达到消元的目的了，但是有一些小小的瑕疵
我们可能碰到一个这样矩阵，有一行全是0,例如这个：
那么我们在步骤1中搜索到主元为0的话，0的任意倍数都是0，会导致第2步无法进行。所以我们需要添加换行的操作，计算方法为：
所以我们把代码逻辑修改成这样：
- 从第一行开始到最后一行
- 找主元：找出第i的主元(第i行第i个数)，若主元为0，把第i行向下换行，直到找到有主元的行。若找不到主元，就开始找下一个
- 消元：消除下面所有行的第i个数（下面每一行减去x倍的第一行来消除第i列）
下面就是高斯消元的主程序：
template <typename T>
bool Matrix<T>::GaussianElimination()
{
Matrix<T> outputMatrix = *this;
/*Gaussian elmiation*/
for(int k=0;k<outputMatrix.m_iRows;k++)
{
/*if all the pivot have been found*/
if(k>=m_iColumns)
{
break;
}
/*exchange rows downward to find the row‘s pivot*/
for(int i=k+1;i<outputMatrix.m_iRows;i++)
{
/*pivot is non-zero*/
if(outputMatrix.m_vecMatrix[k][k] != 0)
{
//T temp = outputMatrix.m_vecMatrix[0][0];
break;
}
else
{
if(i < outputMatrix.m_iRows)
{
outputMatrix.exchangeRows(k,i);
}
}
}
/*if there is no pivot in this row*/
if(outputMatrix.m_vecMatrix[k][k] == 0)
{
break;
}
/*elimination:The rows below pivot row subtract times of pivot row*/
for(int i=k+1;i<outputMatrix.m_iRows;i++)
{
double RowsfirstData = outputMatrix.m_vecMatrix[i][k]/outputMatrix.m_vecMatrix[k][k];//Save the first data of next(k+1) rows
if(RowsfirstData != 0)
{
outputMatrix.m_vecMatrix[i][k]=0;
for(int j=k+1;j<outputMatrix.m_iColumns;j++)
{
outputMatrix.m_vecMatrix[i][j] -= RowsfirstData*outputMatrix.m_vecMatrix[k][j] ;
}
}
}
}
*this = outputMatrix;
return true;
}
高斯-若尔当法
若尔当在高斯消元的基础上加上了回代过程，把矩阵化简成行最简式。我们在高斯消元的基础上加上和回代，方法跟高斯消元相反，用上面的行减下面的行，这里就不详细描述（展开查看代码）
rref()//化简矩阵成行最简
template <typename T>
bool Matrix<T>::rref()
{
Matrix<T> outputMatrix = *this;
int rank=0;//the rank of the matrix, how many columns‘s pivot will it has(-1)
/*Gaussian elmiation*/
for(int k=0;k<outputMatrix.m_iRows;k++)
{
/*if all the pivot elem have been found*/
if(k>=m_iColumns)
{
break;
}
/*exchange rows downward to find the pivot row*/
for(int i=k+1;i<outputMatrix.m_iRows;i++)
{
/*pivot is non-zero*/
if(outputMatrix.m_vecMatrix[k][k] != 0)
{
//T temp = outputMatrix.m_vecMatrix[0][0];
rank++;
break;
}
else
{
if(i < outputMatrix.m_iRows)
{
outputMatrix.exchangeRows(k,i);
}
}
}
/*if there is no pivot in this row*/
if(outputMatrix.m_vecMatrix[k][k] == 0)
{
break;
}
/*elimination:The rows below pivot row subtract times of pivot row*/
for(int i=k+1;i<outputMatrix.m_iRows;i++)
{
double RowsfirstData = outputMatrix.m_vecMatrix[i][k]/outputMatrix.m_vecMatrix[k][k];//Save the first data of next(k+1) rows
if(RowsfirstData != 0)
{
outputMatrix.m_vecMatrix[i][k]=0;
for(int j=k+1;j<outputMatrix.m_iColumns;j++)
{
outputMatrix.m_vecMatrix[i][j] -= RowsfirstData*outputMatrix.m_vecMatrix[k][j] ;
}
}
}
}
/*normalizing:set all pivots to 1*/
for(int i=0;i<outputMatrix.m_iRows;i++)
{
for(int j=0;j<outputMatrix.m_iColumns;j++)
{
if(outputMatrix.m_vecMatrix[i][j] !=0 )//pivot has been foound
{
double pivot = outputMatrix.m_vecMatrix[i][j];//get pivot
for(int k=i;k<outputMatrix.m_iColumns;k++)
{
outputMatrix.m_vecMatrix[i][k] /=pivot;
}
break;
}
}
}
/*Back substitution*/
for(int i = rank;i>=1;i--)
{
/*find a first non-zero elem (It is pivot)*/
for(int j=0;j<outputMatrix.m_iColumns;j++)
{
double times=0;
if(outputMatrix.m_vecMatrix[i][j] !=0)//pivot found
{
for(int l=i-1;l>=0;l--)
{
times = outputMatrix.m_vecMatrix[l][j]/outputMatrix.m_vecMatrix[i][j];
for(int k=j;k<outputMatrix.m_iColumns;k++)//tims of this row subtract by each columns in upon row
{
outputMatrix.m_vecMatrix[l][k] -= times*outputMatrix.m_vecMatrix[i][k];
}
}
break;
}
}
}
*this = outputMatrix;
return true;
}
rrefmovie()//化简矩阵成行最简，并打印过程
template <typename T>
bool Matrix<T>::rrefmovie()
{
Matrix<T> outputMatrix = *this;
int rank=0;//the rank of the matrix, how many columns‘s pivot will it has(-1)
/*Gauss elmiation*/
cout<<"Gauss elimination:"<<endl;
outputMatrix.printfAll();
for(int k=0;k<outputMatrix.m_iRows;k++)
{
/*If all the pivot elem have been found*/
if(k>=m_iColumns)
{
break;
}
/*Exchange rows downward to find the pivot row*/
for(int i=k+1;i<outputMatrix.m_iRows;i++)
{
/*Pivot is non-zero*/
if(outputMatrix.m_vecMatrix[k][k] != 0)
{
rank++;
break;
}
else
{
if(i < outputMatrix.m_iRows)
{
outputMatrix.exchangeRows(k,i);
}
}
if(k!=i)
{
cout<<"row"<<k+1<<" exchange row"<<i+1<<endl;//Debug
outputMatrix.printfAll();
}
}
/*If there is no pivot in this row*/
if(outputMatrix.m_vecMatrix[k][k] == 0)
{
break;
}
/*Elimination:The rows below pivot row subtract times of pivot row*/
for(int i=k+1;i<outputMatrix.m_iRows;i++)
{
double RowsfirstData = outputMatrix.m_vecMatrix[i][k]/outputMatrix.m_vecMatrix[k][k];//Save the first data of next(k+1) rows
if(RowsfirstData != 0)
{
outputMatrix.m_vecMatrix[i][k]=0;
for(int j=k+1;j<outputMatrix.m_iColumns;j++)
{
outputMatrix.m_vecMatrix[i][j] -= RowsfirstData*outputMatrix.m_vecMatrix[k][j] ;
}
}
cout<<"row"<<i+1<<" - "<<RowsfirstData<<"*"<<"row"<<k+1<<endl;//Debug
outputMatrix.printfAll();
}
}
/*Normalizing:set all rows pivot to 1*/
for(int i=0;i<outputMatrix.m_iRows;i++)
{
for(int j=0;j<outputMatrix.m_iColumns;j++)
{
if(outputMatrix.m_vecMatrix[i][j] !=0 )//pivot has been foound
{
double pivot = outputMatrix.m_vecMatrix[i][j];//get pivot
for(int k=i;k<outputMatrix.m_iColumns;k++)
{
outputMatrix.m_vecMatrix[i][k] /=pivot;
}
cout<<"row"<<i+1<<" / "<<pivot<<endl;//Debug
outputMatrix.printfAll();//Debug
break;
}
}
}
/*Back substitution*/
cout<<"Back substitution:"<<endl;
for(int i = rank;i>=1;i--)
{
/*find a first non-zero elem (It is pivot)*/
for(int j=0;j<outputMatrix.m_iColumns;j++)
{
double times=0;
if(outputMatrix.m_vecMatrix[i][j] !=0)//pivot found
{
for(int l=i-1;l>=0;l--)
{
times = outputMatrix.m_vecMatrix[l][j]/outputMatrix.m_vecMatrix[i][j];
for(int k=j;k<outputMatrix.m_iColumns;k++)//tims of this row subtract by each columns in upon row
{
outputMatrix.m_vecMatrix[l][k] -= times*outputMatrix.m_vecMatrix[i][k];
}
cout<<"row"<<l+1<<" - "<<times<<"*"<<"row"<<i+1<<endl;
outputMatrix.printfAll();
}
break;
}
}
}
*this = outputMatrix;
return true;
}
使用我们开始的矩阵测试：
Matrix<double> matrix(3,3);
matrix.setSpecifiedElem(0,0,1);
matrix.setSpecifiedElem(0,1,2);
matrix.setSpecifiedElem(0,2,3);
matrix.setSpecifiedElem(1,0,2);
matrix.setSpecifiedElem(1,1,2);
matrix.setSpecifiedElem(1,2,2);
matrix.setSpecifiedElem(2,0,4);
matrix.setSpecifiedElem(2,1,5);
matrix.setSpecifiedElem(2,2,6);
matrix.printfAll();
matrix.rrefmovie();
matrix.printfAll();
system("pause");
结果：
时间： 2024-10-12 22:30:24