目录
一、 引言…………………………………… 3
二、 算法描述……………………………… 3
三、 试验及分析…………………………… 5
四、 自己的特色…………………………… 7
五、 总结…………………………………… 13
六、 程序清单……………………………… 13
一、引言
在我们的实际生活和生成当中会出现大量的方程组，在这些方程组需要求解，如果是要用人工的去求解的时候无疑是一个庞大的工作，而且很能避免错误的出现，这样会浪费大量的人力物力和时间。于是我们今天在前辈们关于方程组的求解总结出来的高斯消去法和直接三角分解法（LU分解）的基础之上借助编程过程中的循环来让计算机对方程组经行求解。这样可以节省大量的时间同时避免了错误的出现。下面将是我们本着这个目的去实现高斯消去法和直接三角分解法（LU分解）的过程。虽然有很多的难题还没有解决，但是这是我们小组共同努力的成果。
二、算法描述
1、高斯列主元消去法的原理
Gauss消去法的基本思想是一次用前面的方程消去后面的未知数，从而将方程组化为等价形式：
这个过程就是消元，然后再回代就好了。具体过程如下：
对于 ，若 依次计算
然后将其回代得到：
以上是高斯消去。
但是高斯消去法在消元的过程中有可能会出现 的情况，这时消元就无法进行了，即使主元数 但是很小时，其做除数，也会导致其他元素数量级的严重增长和舍入误差的扩散。因此，为了减少误差，每次消元选取系数矩阵的某列中绝对值最大的元素作为主元素。然后换行使之变到主元位置上，再进行销元计算。即高斯列主元消去法。
2、直接三角分解法（LU分解）的原理
先将矩阵A直接分解为 则求解方程组的问题就等价于求解两个三角形方程组。
直接利用矩阵乘法，得到矩阵的三角分解计算公式为：
由上面的式子得到矩阵A的LU分解后，求解Ux=y的计算公式为
以上为LU分解法。
三、实验及分析
高斯列组消元法：
测试数据1：3 4
1 2 3 6
2 2 3 7
1 4 4 9
答案：1 1 1
测试数据2：3 4
0 0 1 3
3 0 0 2
0 2 0 1
答案：0.667 0.5 3
LU算法：
实例一：
实例二：
、
实例三：
四、自己的特色
1用高斯算法解一般的方程组
实例一：
3 5
10 200 300 4 233
22 22 44 1 22
1 20 30 0.4 233
实例二：
4 6
100 34 345 334 11 552
12 0.2 0.12 12 23 344
122 53 12 43 23 45
123 234 0.3 34 345 34
实例三：
3 4
1 22 676 345
12 24 124 466
1 22 676 100
实例四：
4 5
21 345 67 15 456
0.12 35 02.4 234 3455
122 34 13 345 234
12 12 0.1 24 24
实例五：
5 6
1 2 3 4 6 7
1 2 5 12 4 56
12 2 5 11 2 122
32 11 2 4 6 23
2 4 8 16 10 63
实例六：
4 3
1 3 45
32 12 24
42 0.1 23
33 15 60
实例七：
3 3
1 2 4
1 3 5
2 5 9
实例八：
4 4
1 4 7 3
2 8 14 6
3 12 21 9
3 3 6 2
1 两种高斯消去法比较
(1)列主元高斯消去法
(2)全主元高斯消去法
五、总结
通过这次实习我们学到了很多同样也暴露出了很多的问题：
1、这从的方程组求解的两种方法不像前一次迭代的那样好实现，所以在很多的方面大家写的代码又点不完善的地方，有几个同学本来的基础就是有点弱，所以全部的程序没有写完，有的写完了但是会出现些错误，这样暴露出来在我们编程能力不一样，由于上课大家在一起讲解的时候时间有限，这样我们思考要是还是一个人单独作战的时候是不是会使大家的水平更加的层次不起。所以我们决定进行小组内再分组的方案。有两个人组成一组，形成一个变成能力比较好的人带一个比较差的人，这样大家在课下编程的时候要一起行动，相互的帮助可能比我们在课上汇报讨论的时候更能解决实质问题。
2、由于这次大家都对原理比较清楚，在私下也有点交流，所以在两种方法的主要实现过程中没有太大的问题大大加速了汇报的进度。在最后的时候我们进行的问题的集中解决，并且在最后的时候请别的小组的实现的比较好的同学给我们指导，并且相互交流。明白别的小组对我们小组出现的问题是怎么解决的。然后在最后得出了一些大家都有的一些共有的缺陷，进行的商讨。
3、在即上一次的分工之后，我们这次决定相互的交换一下工作，然后让一些编程能力相对比较弱的同学参与到了讨论后的代码修改当中。我们最终的程序或许没有解决我们最初讨论出的问题，但是我想大家参与了就是一种最大的收获。这次的分工是：由彭亚琼、刘宇、方正、陈俊桦、韩学武五位同学负责修改程序，王艳琴和向胜男负责细化许多上台汇报的工作，张钟文和杨耀鹏写word文档，肖大军和梅旭负责制作PPT。
六、程序清单
高斯列主消元
#include<stdio.h>
#include<string.h>
#include<math.h>
#define e 0.00000001
const int M=1000;
const int N=1000;
double a[N][M];
double b[M];
double pt(int n,int m)//矩阵的输出
{
for(int i=0;i<n;i++)
{
for(int j=0;j<m;j++)
printf("%lf ",a[i][j]);
printf("\n");
}
}
int jihuan(double a[N][M],int n,int m,int i)
{
int mi=i;//最大行首
//printf("%lf\n",a[mi][mi]);
int x;
for(int j=i+1;j<n;j++)
if(fabs(a[j][i])>fabs(a[mi][i]))mi=j;//列不动
//printf("##%d\n",mi);
//printf("%.6f\n",a[mi][i]);
if(mi!=i)
for(int j=0;j<m;j++)//换行
{
x=a[mi][j];
a[mi][j]=a[i][j];
a[i][j]=x;
}
//printf("##\n");
//pt(n,m);
//printf("##\n");
return 0;
}
int main()
{
int n,m,zhi;
double t;
printf("请输入增广矩阵N,M\n");
while(scanf("%d %d",&n,&m)!=EOF)
{
for(int i=0;i<n;i++)
for(int j=0;j<m;j++)
scanf("%lf",&a[i][j]);
//pt(n,m);
for(int i=0;i<n-1;i++)
{
jihuan(a,n,m,i);//找到不为0的因子,交换行首最大
for(int j=i+1;j<n;j++)
{
if(fabs(a[j][i])>e)
{
t=a[j][i]/a[i][i];//对a[i][i]处理
for(int k=i;k<m;k++)
{
a[j][k]=a[j][k]-a[i][k]*t;
}
}
}
pt(n,m);printf("\n");
}
int j;
for(j=0;j<m;j++)
if(fabs(a[n-1][j])>e)break;//判断最后行是否全为0
//printf("%lf %lf",a[n-1][m-2],a[n-1][m-1]);
if(fabs(a[n-1][m-2])<e&&fabs(a[n-1][m-1])>e){printf("方程无解\n");continue;}
//printf("%d %d",j,m);
if(j==m)
{
for(int i=n-1;i>-1;i--)//其他行
if(fabs(a[i][m-2])<e)n--;
}
if(n<m-1)//判断矩阵的秩
{
printf("方程有无穷多解\n");
printf("化简后的矩阵为：\n");
pt(n,m);
printf("请输入拓展矩阵N,M\n");
continue;
}
else if(n>m-1)
{
printf("方程无解\n");
printf("化简后的矩阵为：\n");
pt(n,m);
printf("请输入拓展矩阵N,M\n");
continue;
}
for(int i=n-1;i>0;i--)
for(int j=i-1;j>-1;j--)
{
t=a[j][i]/a[i][i];
//printf("%.2f\n",t);
a[j][i]=a[j][i]-a[i][i]*t;//对a[i][i]处理
a[j][m-1]-=a[i][m-1]*t;
}
pt(n,m);
printf("\n");
printf("方程的解为：");
for(int i=0;i<n;i++)
printf("%.3lf ",a[i][m-1]/a[i][i]);
printf("\n请输入拓展矩阵N,M\n");
}
return 0;
}
高斯全主消元
#include<iostream>
#include<string.h>
#include<cstdio>
#include<cmath>
#include<iomanip>
using namespace std;
#define exp 0.01
double linear[101][101];
double pointer[101];//////////////////////////////////存放每行下标
double x[101];
double root[100];
int n;
double iteration_root(int i)
//求方程组中第i个未知数的解
{
if(i==n)
return linear[n][n+1];
else
{
double sum=linear[i][n+1];
for(int j=i;j<n;j++)
sum-=linear[i][j+1]*iteration_root(j+1);
return sum;
}
}
void Gauss_elimination()
//Gauss_elimination化为上三角矩阵
{
int t,i,j,k;
int tag=0;
for(k=1;k<=n;k++)
{
int max_index1=k,max_index2=k;
for(i=k;i<=n;i++)
{
for(int j=k;j<=n;j++)
if(fabs(linear[i][j])>fabs(linear[max_index1][max_index2]))
{
max_index1=i;
max_index2=j;
}
}
//if(fabs(linear[max_index1][max_index2])>=exp)
for(i=1;i<=n+1;i++)
{
pointer[i]=linear[max_index1][i];
linear[max_index1][i]=linear[k][i];
linear[k][i]=pointer[i];
}
for(i=1;i<=n+1;i++)
{
pointer[i]=linear[i][max_index2];
linear[i][max_index2]=linear[i][k];
linear[i][k]=pointer[i];
}
for(i=k;i<=n;i++)
{
if(i==k)
{
double w=linear[i][i];
for(j=k;j<=n+1;j++)
linear[i][j]/=w;
}
else
{
double w=linear[i][k];
for(j=k;j<=n+1;j++)
linear[i][j]=linear[i][j]-w*linear[k][j];
}
}
tag++;
cout<<"经过"<<tag<<"次变换线性方程组对应的线性矩阵为:"<<endl;
for(int ii=1;ii<=n;ii++)
{
for( t=1;t<=n;t++)
cout<<setw(15)<<linear[ii][t];
cout<<"="<<setw(15)<<linear[ii][t]<<endl;
}
getchar();
}
}
int main()
{
char varible[100];
for(int i=0;i<100;i++)
{
varible[i]='x';
}
int i,j;
cout<<"请输入系数矩阵的阶数n"<<endl;
while(cin>>n)
{
cout<<"请输入系数矩阵"<<endl;
for(i=1;i<=n;i++)
{
for(j=1;j<=n+1;j++)
{
scanf("%lf",&linear[i][j]);
}
}
for(i=1;i<n+1;i++)
{
linear[n+1][i]=i;
}
getchar();
cout<<"输入的线性方程组为:"<<endl;
for(i=1;i<=n;i++)
{
for(j=1;j<=n;j++)
cout<<linear[i][j]<<"*"<<varible[j-1]<<j<<"+";
cout<<"\b="<<linear[i][j]<<endl;
}
Gauss_elimination();
for(i=1;i<=n;i++)
root[i]=iteration_root(i);
for(i=1;i<=n;i++)
{
int t=linear[n+1][i];
x[t]=root[i];
}
cout<<"线性方程组的解为:"<<endl;
for(i=1;i<=n;i++)
cout<<varible[i-1]<<i<<"="<<setw(10)<<setprecision(8)<<x[i]<<endl;
cout<<endl<<endl<<"请输入系数矩阵的阶数n"<<endl;
}
return 0;
}
LU算法：
#include<iostream>
#include<string.h>
using namespace std;
double m[20][20],lu[20][20];//系数，LU
double root[20],x[20],y[20];
int n;
void LU()//化解
{
int i,j,k,t;
double num;
memset(lu,0,sizeof(lu));
for(i=0;i<n;i++)//U的第一行，L的第一列处理
{
lu[0][i]=m[0][i];
if(i>0)
lu[i][0]=m[i][0]/lu[0][0];
}
for(i=1;i<n;i++)
{
for(j=i;j<n;j++)//U的第i行
{
num=0;
for(t=0;t<i;t++)
{
//if(i-1==t) {num+=lu[t][j];continue;}
num+=lu[i][t]*lu[t][j];
}
//if(i==1) cout<<num<<endl;
lu[i][j]=m[i][j]-num;
}
for(j=i+1;j<n;j++)//L的第i列
{
num=0;
for(k=0;k<i;k++)
num+=lu[j][k]*lu[k][i];
lu[j][i]=(m[j][i]-num)/lu[i][i];
}
}
cout<<endl;
cout<<"分解后的LU矩阵为："<<endl;
for(i=0;i<n;i++)
{
for(j=0;j<n;j++)
cout<<lu[i][j]<<" ";
cout<<endl;
}
cout<<endl;
}
double iteration(int i)
//求方程组中第i个未知数的解（递推）
{
double sum;
if(i==0)
return root[0];
else
{
double sum=root[i];
for(int j=i-1;j>=0;j--)
{
sum-=lu[i][j]*iteration(j);
//cout<<lu[i][j]<<" ";
}
return sum;
}
}
double iteration_root(int i)
//求方程组中第i个未知数的解
{
double sum;
if(i==n-1)
return y[n-1]/lu[i][i];
else
{
double sum=y[i];
for(int j=i+1;j<n;j++)
{
sum=(sum-lu[i][j]*iteration_root(j));
}
sum=sum/lu[i][i];
return sum;
}
}
int main()
{
int i,j,k,t;
double s;
memset(m,0,sizeof(m));
memset(root,0,sizeof(root));
memset(x,0,sizeof(x));
memset(y,0,sizeof(y));
cout<<"请输入系数矩阵的阶数"<<endl;
cin>>n;
cout<<"请输入系数矩阵"<<endl;
for(i=0;i<n;i++)
for(j=0;j<n;j++)
cin>>m[i][j];
cout<<"请输入结果数组"<<endl;
for(i=0;i<n;i++)
cin>>root[i];
for(i=0;i<n;i++)//如果系数矩阵对角线上存在为0的元素，则对矩阵进行旋转
{
s=0;
if(m[i][i]==0)
{
for(j=0;j<n;j++)
{
if(m[j][i]>s)
{
s=m[j][i];
t=j;
}
}
for(k=0;k<n;k++)
{
s=m[t][k];
m[t][k]=m[i][k];
m[i][k]=s;
}
s=root[t];
root[t]=root[i];
root[i]=s;
}
}
LU();
for(j=0;j<n;j++)
{
y[j]=iteration(j);
cout<<"Y"<<j<<"="<<y[j]<<endl;
}
cout<<endl;
for(k=n-1;k>=0;k--)
x[k]=iteration_root(k);
cout<<"方程解为:"<<endl;
for(k=0;k<n;k++)
cout<<"X"<<k<<"="<<x[k]<<endl;
cout<<endl;
}
#include<iostream> #include<string.h> using namespace std; double m[20][20],lu[20][20];//系数，LU double root[20],x[20],y[20]; int n; void LU()//化解 { int i,j,k,t; double num; memset(lu,0,sizeof(lu)); for(i=0;i<n;i++)//U的第一行，L的第一列处理 { lu[0][i]=m[0][i]; if(i>0) lu[i][0]=m[i][0]/lu[0][0]; } for(i=1;i<n;i++) { for(j=i;j<n;j++)//U的第i行 { num=0; for(t=0;t<i;t++) { //if(i-1==t) {num+=lu[t][j];continue;} num+=lu[i][t]*lu[t][j]; } //if(i==1) cout<<num<<endl; lu[i][j]=m[i][j]-num; } for(j=i+1;j<n;j++)//L的第i列 { num=0; for(k=0;k<i;k++) num+=lu[j][k]*lu[k][i]; lu[j][i]=(m[j][i]-num)/lu[i][i]; } } cout<<endl; cout<<"分解后的LU矩阵为："<<endl; for(i=0;i<n;i++) { for(j=0;j<n;j++) cout<<lu[i][j]<<" "; cout<<endl; } cout<<endl; } double iteration(int i) //求方程组中第i个未知数的解（递推） { double sum; if(i==0) return root[0]; else { double sum=root[i]; for(int j=i-1;j>=0;j--) { sum-=lu[i][j]*iteration(j); //cout<<lu[i][j]<<" "; } return sum; } } double iteration_root(int i) //求方程组中第i个未知数的解 { double sum; if(i==n-1) return y[n-1]/lu[i][i]; else { double sum=y[i]; for(int j=i+1;j<n;j++) { sum=(sum-lu[i][j]*iteration_root(j)); } sum=sum/lu[i][i]; return sum; } } int main() { int i,j,k,t; double s; memset(m,0,sizeof(m)); memset(root,0,sizeof(root)); memset(x,0,sizeof(x)); memset(y,0,sizeof(y)); cout<<"请输入系数矩阵的阶数"<<endl; cin>>n; cout<<"请输入系数矩阵"<<endl; for(i=0;i<n;i++) for(j=0;j<n;j++) cin>>m[i][j]; cout<<"请输入结果数组"<<endl; for(i=0;i<n;i++) cin>>root[i]; for(i=0;i<n;i++)//如果系数矩阵对角线上存在为0的元素，则对矩阵进行旋转 { s=0; if(m[i][i]==0) { for(j=0;j<n;j++) { if(m[j][i]>s) { s=m[j][i]; t=j; } } for(k=0;k<n;k++) { s=m[t][k]; m[t][k]=m[i][k]; m[i][k]=s; } s=root[t]; root[t]=root[i]; root[i]=s; } } LU(); for(j=0;j<n;j++) { y[j]=iteration(j); cout<<"Y"<<j<<"="<<y[j]<<endl; } cout<<endl; for(k=n-1;k>=0;k--) x[k]=iteration_root(k); cout<<"方程解为:"<<endl; for(k=0;k<n;k++) cout<<"X"<<k<<"="<<x[k]<<endl; cout<<endl; }
#include<stdio.h> #include<string.h> #include<math.h> #define e 0.00000001 const int M=1000; const int N=1000; double a[N][M]; double b[M]; double pt(int n,int m)//矩阵的输出 { for(int i=0;i<n;i++) { for(int j=0;j<m;j++) printf("%lf ",a[i][j]); printf("\n"); } } int jihuan(double a[N][M],int n,int m,int i) { int mi=i;//行首 //printf("%lf\n",a[mi][mi]); int x; for(int j=i+1;j<n;j++) if(fabs(a[j][i])>fabs(a[mi][i]))mi=j;//列不动 //printf("##%d\n",mi); //printf("%.6f\n",a[mi][i]); if(mi!=i) for(int j=0;j<m;j++)//换行 { x=a[mi][j]; a[mi][j]=a[i][j]; a[i][j]=x; } //printf("##\n"); //pt(n,m); //printf("##\n"); return 0; } int main() { int n,m,zhi; double t; printf("请输入增广矩阵N,M\n"); while(scanf("%d %d",&n,&m)!=EOF) { for(int i=0;i<n;i++) for(int j=0;j<m;j++) scanf("%lf",&a[i][j]); //pt(n,m); for(int i=0;i<n-1;i++) { jihuan(a,n,m,i);//找到不为0的因子,交换行首最大 for(int j=i+1;j<n;j++) { if(fabs(a[j][i])>e) { t=a[j][i]/a[i][i];//对a[i][i]处理 for(int k=i;k<m;k++) a[j][k]=a[j][k]-a[i][k]*t; } } pt(n,m);printf("\n"); } int j; for(j=0;j<m;j++) if(fabs(a[n-1][j])>e)break;//判断最后行是否全为0 //printf("%lf %lf",a[n-1][m-2],a[n-1][m-1]); //for(j=n-1;j>-1;j--) if(fabs(a[j][m-2])<e&&fabs(a[j][m-1])>e){printf("方程无解\n");continue;} //printf("%d %d",j,m); //if(ki==1)continue; if(j==m) { for(int i=n-1;i>-1;i--)//其他行 if(fabs(a[i][m-2])<e)n--; } if(n<m-1)//判断矩阵的秩 { printf("方程有无穷多解\n"); printf("化简后的矩阵为：\n"); pt(n,m); printf("请输入拓展矩阵N,M\n"); continue; } else if(n>m-1) { printf("方程无解\n"); printf("化简后的矩阵为：\n"); pt(n,m); printf("请输入拓展矩阵N,M\n"); continue; } for(int i=n-1;i>0;i--) for(int j=i-1;j>-1;j--) { t=a[j][i]/a[i][i]; //printf("%.2f\n",t); a[j][i]=a[j][i]-a[i][i]*t;//对a[i][i]处理 a[j][m-1]-=a[i][m-1]*t; } pt(n,m); printf("\n"); printf("方程的解为："); for(int i=0;i<n;i++) printf("%.3lf ",a[i][m-1]/a[i][i]); printf("\n请输入拓展矩阵N,M\n"); } return 0; }
#include<iostream> #include<string.h> #include<cstdio> #include<cmath> #include<iomanip> using namespace std; #define exp 0.01 double linear[101][101]; double pointer[101];//////////////////////////////////存放每行下标 double x[101]; double root[100]; int n; double iteration_root(int i) //求方程组中第i个未知数的解 { if(i==n) return linear[n][n+1]; else { double sum=linear[i][n+1]; for(int j=i;j<n;j++) sum-=linear[i][j+1]*iteration_root(j+1); return sum; } } void Gauss_elimination() //Gauss_elimination化为上三角矩阵 { int t,i,j,k; int tag=0; for(k=1;k<=n;k++) { int max_index1=k,max_index2=k; for(i=k;i<=n;i++) { for(int j=k;j<=n;j++) if(fabs(linear[i][j])>fabs(linear[max_index1][max_index2])) { max_index1=i; max_index2=j; } } //if(fabs(linear[max_index1][max_index2])>=exp) for(i=1;i<=n+1;i++) { pointer[i]=linear[max_index1][i]; linear[max_index1][i]=linear[k][i]; linear[k][i]=pointer[i]; } for(i=1;i<=n+1;i++) { pointer[i]=linear[i][max_index2]; linear[i][max_index2]=linear[i][k]; linear[i][k]=pointer[i]; } /*else if(linear[max_index1][max_index2]!=0) { for(i=k;i<=m+1;i++) linear[i][k]*=(10*exp/linear[k][k]); }*/ //double w=linear[k][k]; for(i=k;i<=n;i++) { if(i==k) { double w=linear[i][i]; for(j=k;j<=n+1;j++) linear[i][j]/=w; } else { double w=linear[i][k]; for(j=k;j<=n+1;j++) linear[i][j]=linear[i][j]-w*linear[k][j]; } } tag++; cout<<"经过"<<tag<<"次变换线性方程组对应的线性矩阵为:"<<endl; for(int ii=1;ii<=n;ii++) { for( t=1;t<=n;t++) cout<<setw(15)<<linear[ii][t]; cout<<"="<<setw(15)<<linear[ii][t]<<endl; } getchar(); } } int main() { char varible[100]; for(int i=0;i<100;i++) { varible[i]='x'; } int i,j; cout<<"请输入系数矩阵的阶数n"<<endl; while(cin>>n) { cout<<"请输入系数矩阵"<<endl; for(i=1;i<=n;i++) { for(j=1;j<=n+1;j++) { scanf("%lf",&linear[i][j]); } } for(i=1;i<n+1;i++) { linear[n+1][i]=i; } getchar(); cout<<"输入的线性方程组为:"<<endl; for(i=1;i<=n;i++) { for(j=1;j<=n;j++) cout<<linear[i][j]<<"*"<<varible[j-1]<<j<<"+"; cout<<"\b="<<linear[i][j]<<endl; } Gauss_elimination(); for(i=1;i<=n;i++) root[i]=iteration_root(i); for(i=1;i<=n;i++) { int t=linear[n+1][i]; x[t]=root[i]; } cout<<"线性方程组的解为:"<<endl; for(i=1;i<=n;i++) cout<<varible[i-1]<<i<<"="<<setw(10)<<setprecision(8)<<x[i]<<endl; cout<<endl<<endl<<"请输入系数矩阵的阶数n"<<endl; } return 0; }