数值计算（一）之解线性方程组（高斯消去法，列选主元消去法，全选主元消去法，杜立特尔分解，克洛特分解，乔里斯基分解）
解线性方程组即解一个多元一次方程组，例如
消去法
原理
没有学过高级的解法也没关系，凭借我们初高中的知识足以解决这个问题
这是一个多元一次方程组，拥有n个未知量，也有n方程
我们可以通过对一个方程组进行同等变形后与其他方程组相加减，消去其他未知量
最后变为一个仅含有一个未知量的方程，求解后逐次回带，接触其他未知数
当我们接触了线性代数后，我们可以把这n个方程组的系数提取出来，与等号右边的值共同组成一个增广矩阵，对这个矩阵进行初等行列变换，然后求解。实质上和我们初高中的消元法是一样的。
例如
我们构造这样一个 
数值计算这里把解线性方程组的方法系统化，分为消去法，分解法，迭代法。
1分解法
分解法是将我们上面构造的矩阵，进行初等行列变换后，变成单位上三角矩阵或者单位下三角矩阵，然后求解出一个未知数后逐个回带。
高斯消去法就是把矩阵变换为单位上三角矩阵
求解出X3以后逐个回带
列选消去法是在上面的基础上进行改进，每次选择对角线下那一列元素绝对值最大的元素，调整到对角线部分，在进行处理，用以减小误差。
全选消去法又是在列选消去法的基础上改进，每次选择从对角线元素开始的右下角矩形部分的最大值，把它放在对角线上，用以减小误差
由于全选消去法可能会交换列，进行列变换，改变变量的位置，所以用计算机求解时需要记录交换的情况，最后把未知数恢复为X1-Xn的顺序
高斯消去法
//高斯消去法
void GusElimination(double *p,int row)
{
double coef;//系数，用于归一化和行处理
int col=row+1;//保存增广矩阵的列数
int i,j,k;//循环变量，指示行和列
for(k=0;k<row;k++)//有多少个未知量即循环多少次
{
coef=*(p+k*col+k);
for(j=0;j<col;j++)//对增广矩阵的第k行做归一化处理
*(p+k*col+j)/=coef;
for(i=k+1;i<row;i++)//消去k+1-row行对角线以下的元素
{
coef=(-1)*(*(p+i*col+k))/(*(p+k*col+k));//求与第k行的系数
for(j=0;j<col;j++)
*(p+i*col+j)+=(coef*(*(p+k*col+j)));//对k+1-row行所有元素都加上coef倍的第k行
}//End for-i
}//End for-k
//求得单位上三角矩阵后，从下往上求解未知量
for(i=row-1;i>=0;i--)
for(j=row-1;j>i;j--)
*(p+i*col+row)-=( *(p+j*col+row) ) * ( *(p+i*col+j) );
//求解结果保存在增广矩阵的最后一列
} 列选消去法
//列选主元消去法
void ColumnElimination(double *p,int row)
{
double coef,tmp,Max;
int Maxrow,col=row+1;
int i,j,k;
for(k=0;k<row;k++)
{
Max=fabs(*(p+k*col+k));
Maxrow=k;
for(i=k+1;i<row;i++)//列选过程，找到当列绝对值最大的元素
if( fabs(*(p+i*col+k)) > Max){
Max=*(p+i*col+k); Maxrow=i;
}//End if
for(j=k;j<col;j++)
{
tmp = *(p+k*col+j);
*(p+k*col+j)=*(p+Maxrow*col+j);
*(p+Maxrow*col+j)=tmp;
}//End for-j
coef=*(p+k*col+k);//归一化处理，可以和上面步骤合并
for(j=k;j<col;j++)
*(p+k*col+j)/=coef;
for(i=k+1;i<row;i++)//消去对角线下的元素
{
coef= (-1) * ( *(p+i*col+k) ) / ( *(p+k*col+k) );
for(j=k;j<col;j++)
*(p+i*col+j) += coef * ( *(p+k*col+j) );
}//End for-i
}//End for-k
//求得单位上三角矩阵后，从下往上求解未知量
for(i=row-1;i>=0;i--)
for(j=row-1;j>i;j--)
*(p+i*col+row)-=( *(p+j*col+row) ) * ( *(p+i*col+j) );
//求解结果保存在增广矩阵最后一列
} 全选主元消去法
//全选主元消去法
void AllElimination(double *p,int row)
{
double coef,tmp,Max;//coef表系数，tmp临时变量，Max用来保存绝对值最大元素
int MaxCol,MaxRow,col=row+1;//MaxCol和MaxRow用来保存绝对值最大元素所在的行和列
int i,j,k;
int *a=(int *)malloc(row * sizeof(int));//申请一个临时数组用来保存交换记录
for(i=0;i<row;i++)//对a数组进行初始化
*(a+i)=i;
for(k=0;k<row;k++)
{
Max=fabs(*(p+k*col+k));
coef=*(p+k*col+k);
MaxCol=k; MaxRow=k;
for(i=k;i<row;i++)//全选的过程
{
for(j=k;j<row;j++)
{
if( fabs(*(p+i*col+j)) >Max){
Max = fabs(*(p+i*col+j));
coef=*(p+i*col+j);
MaxRow=i; MaxCol=j;
}// End if
}//End for-j
}//End for-i
if(MaxCol != k)//交换k和MaxCol两列，行i的范围从0到Row-1
{
for(i=0;i<row;i++)
{
tmp = *(p+i*col+k);
*(p+i*col+k)=*(p+i*col+MaxCol);
*(p+i*col+MaxCol)=tmp;
}//End for-i
//交换两列时需要记录
tmp = *(a+k);//记录交换的是哪两列
*(a+k)=*(a+MaxCol);
*(a+MaxCol)=tmp;
}//End if
for(j=k;j<col;j++)
{//一起处理需要先进行列交换(如果需要)
if(MaxRow != k)
{
tmp=*(p+k*col+j);
*(p+k*col+j)=*(p+MaxRow*col+j);
*(p+MaxRow*col+j)=tmp;
}//End if
*(p+k*col+j)/=coef;
}//End for-j
for(i=k+1;i<row;i++)//消去对角线下的元素
{
coef= (-1) * ( *(p+i*col+k) ) / ( *(p+k*col+k) );
for(j=k;j<col;j++)
*(p+i*col+j) += coef * ( *(p+k*col+j) );
}
}
//求得单位上三角矩阵后，从下往上求解未知量
for(i=row-1;i>=0;i--)
for(j=row-1;j>i;j--)
*(p+i*col+row)-=( *(p+j*col+row) ) * ( *(p+i*col+j) );
for(i=0;i<row-1;i++)
{//对a进行冒泡排序，同时求解结果列进行相同操作
for(j=0;j<row-i-1;j++)//通过排序来复原未知量原来的顺序
{
if( *(a+j) > *(a+j+1) )
{
tmp = *(a+j);
*(a+j)=*(a+j+1);
*(a+j+1)=tmp;
tmp = *(p+j*col+row);
*(p+j*col+row)= *(p+(j+1)*col+row);
*(p+(j+1)*col+row)=tmp;
}//End if
}//End for-j
}//End for-i
free(a);
} 分解法
原理
消去法是增广矩阵化为单位上三角矩阵求解，除此之外，还有另外一种解题思路
那就是：将矩阵A分解成两个三角形矩阵L和U的乘积，其中L为下三角矩阵，U为上三角矩阵，则
可以证明，当系数矩阵的各阶顺序主子式不为0时，这样的分为有且仅有一组
在这种思路下，就有三种分解方法，
Doolitte分解法：L是单位下三角矩阵，U是上三角矩阵
Crout分解法：L是下三角矩阵，U是单位上三角矩阵
Cholesky分解法：L和U互为转置矩阵
要进行这样的分解，系数矩阵A就要满足是对称正定矩阵
对称：矩阵中的元素关于对角线对称
正定：矩阵所有的特征值都大于0
三角分解法不需要特别记住什么公式，只需要知道矩阵乘法的运算方式
进行一次逆运算，先计算某一个矩阵的行，在计算另外一个矩阵的列，以此往复。
杜立特尔分解法
//杜立特尔分解法,分解为单位下三角和上三角矩阵
void Doolittle(double *p,int row)
{
int i,j,k,col=row+1;//增广矩阵的列
double sum;//保存求和结果
for(k=0;k<row;k++)//K表示循环次数，也指示矩阵的行和列，即沿对角线移动
{//先计算U矩阵的K行
for(j=k;j<row;j++)//控制U矩阵列的移动
{//从对角线后面开始求,所以下标从k开始
sum=0;//L矩阵固定K行
for(i=0;i<k;i++)//乘到对角线前一个数字
sum+= ( *(p+k*col+i) ) * ( *(p+i*col+j) );
*(p+k*col+j) = ( *(p+k*col+j) ) - sum;
}//End for-j
//计算L矩阵的列
for(i=k+1;i<row;i++)//控制L矩阵行的移动
{//只需要计算对角线下面的元素,所以下标从k+1开始
sum=0;
for(j=0;j<k;j++)//U矩阵固定K列
sum+= ( *(p+i*col+j) ) * ( *(p+j*col+k) );
*(p+i*col+k)=( ( *(p+i*col+k) ) -sum)/( *(p+k*col+k) );
}//End for-i
}//End for-k
for(i=0;i<row;i++) //通过Ly=b求Ux的结果y
{//结果保存在增广矩阵最后一列
sum=0;
for(j=0;j<i;j++)
sum+=( *(p+j*col+row) )*( *(p+i*col+j) );
*(p+i*col+row)-=sum;
}//End for-i
for(i=row-1;i>=0;i--)//通过Ux=y求解x
{
for(j=row-1;j>i;j--)
*(p+i*col+row)-=( *(p+j*col+row) ) * ( *(p+i*col+j) );
( *(p+i*col+row) )/=( *(p+i*col+i) );
}//End for-i
for(i=0;i<row;i++)
{
for(j=0;j<col;j++)
printf("%lf ",*(p+i*col+j));
printf("\n");
}
} 克洛特分解法
//克洛特尔分解法，分解为单位上三角和下三角矩阵
void Crout(double *p,int row)
{
int i,j,k;
int col=row+1;
double sum;
for(k=0;k<row;k++)//先求L矩阵的列
{
for(i=k;i<row;i++)//i控制矩阵L行的移动
{
sum=0;
for(j=0;j<k;j++)//j控制L矩阵列的移动，需要累乘至对角线前一个
sum+= ( *(p+i*col+j) )*( *(p+j*col+k) );
*(p+i*col+k)= ( *(p+i*col+k)-sum );//求L的k列i行
}//End for-i
for(j=k+1;j<row;j++)//求U矩阵的第k行
{//j控制矩阵U列的移动
sum=0;
for(i=0;i<k;i++)//i控制U矩阵行的移动
sum+=( *(p+k*col+i) )*( *(p+i*col+j) );
*(p+k*col+j)=( *(p+k*col+j)-sum )/( *(p+k*col+k) );
}//求U的k行j列
}//End for-k
for(i=0;i<row;i++)
{//通过Ly=b求解结果y
sum=0;
for(j=0;j<i;j++)
sum+=( *(p+i*row+j) )*( *(p+j*col+row) );
*(p+i*col+row)=( (*(p+i*col+row))-sum )/( *(p+i*row+i) );
}//End for-i
for(i=row-1;i>=0;i--)//通过Ux=y求解x
for(j=row-1;j>i;j--)
*(p+i*col+row)-=( *(p+j*col+row) )*( *(p+i*col+j) );
for(i=0;i<row;i++)
{
for(j=0;j<col;j++)
printf("%lf ",*(p+i*col+j));
printf("\n");
}
} 乔里斯基分解法
//乔里斯基分解，分解为两个互为转置矩阵的矩阵
void Cholesky(double *p,int row)
{
int i,j,k;
int col=row+1;
double sum;
for(k=0;k<row;k++)
{
for(j=k;j<row;j++)//先求U的行
{
sum=0;
for(i=0;i<k;i++)
sum+=( *(p+k*col+i)*( *(p+i*col+j) ) );
if(j==k)//对角线部分需要开根号
*(p+k*col+j)=sqrt( ( *(p+k*col+j) )-sum );
else
*(p+k*col+j)=(( *(p+k*col+j))-sum)/( *(p+k*col+k) );
*(p+j*col+k)=*(p+k*col+j);//L和U互为转置，关于对角线对称
}//End for-j
}//End for-k
for(i=0;i<row;i++)
{//通过Ly=b求解结果y
for(j=0;j<i;j++)
*(p+i*col+row)-=( *(p+i*col+j) )*( *(p+j*col+row) );
*(p+i*col+row)/=( *(p+i*col+i) );
}
for(i=row-1;i>=0;i--)
{//通过Ux=y求解结果x
for(j=row-1;j>i;j--)
*(p+i*col+row)-=( *(p+i*col+j) )*( *(p+j*col+row) );
*(p+i*col+row)/=( *(p+i*col+i) );
}
for(i=0;i<row;i++)
{
for(j=0;j<col;j++)
printf("%lf ",*(p+i*col+j));
printf("\n");
}
}
