目录

 

一、 引言…………………………………… 3

 

二、 算法描述……………………………… 3

 

三、 试验及分析…………………………… 4

 

四、 自己的特色…………………………… 6

 

五、 总结…………………………………… 9

 

六、 程序清单……………………………… 9

                                            

 

 

 

 

 

 

一、引言

     在我们的实际生活和生成当中会出现大量的方程组，在这些方程组需要求解，如果是要用人工的去求解的时候无疑是一个庞大的工作，而且很能避免错误的出现，这样会浪费大量的人力物力和时间。于是我们今天在前辈们关于方程组的求解总结出来的高斯消去法和直接三角分解法（LU分解）的基础之上借助编程过程中的循环来让计算机对方程组经行求解。这样可以节省大量的时间同时避免了错误的出现。下面将是我们本着这个目的去实现高斯消去法和直接三角分解法（LU分解）的过程。虽然有很多的难题还没有解决，但是这是我们小组共同努力的成果。

 

二、算法描述

1、高斯列主元消去法的原理

Gauss消去法的基本思想是一次用前面的方程消去后面的未知数，从而将方程组化为等价形式：

 

这个过程就是消元，然后再回代就好了。具体过程如下：

对于 ，若 依次计算

 

然后将其回代得到：

 

以上是高斯消去。

但是高斯消去法在消元的过程中有可能会出现 的情况，这时消元就无法进行了，即使主元数 但是很小时，其做除数，也会导致其他元素数量级的严重增长和舍入误差的扩散。因此，为了减少误差，每次消元选取系数矩阵的某列中绝对值最大的元素作为主元素。然后换行使之变到主元位置上，再进行销元计算。即高斯列主元消去法。

2、直接三角分解法（LU分解）的原理

   先将矩阵A直接分解为 则求解方程组的问题就等价于求解两个三角形方程组。

直接利用矩阵乘法，得到矩阵的三角分解计算公式为：

 

由上面的式子得到矩阵A的LU分解后，求解Ux=y的计算公式为

 

以上为LU分解法。

 

三、实验及分析

高斯

 

3 4

1 2 3 6

2 2 3 7

1 4 4 9

答案：1 1 1

 

3 4

0 0 1 3

3 0 0 2

0 2 0 1

答案：0.667   0.5   3

 

LU算法：

中文书本P103例7:

 

 

 

实例:  系数矩阵0 1

               1 1     结果矩阵 1 1

 

 

 

 

 

四、自己的特色

 

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

 

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
      
       #include
       
        #define zero 1e-5                   //允许的误差范围double x[10000];                  //保存每一次的计算值double fr(double x){    return x-exp(-x);             //测试函数f(x)}double fd(double x){    return exp(-x);               //迭代法迭代式}double fn(double x)              //牛顿法迭代公式{    return x-(x-exp(-x))/(1+x);}int erfen(double a,double b){    int i=0,j;    double x1,x2,x3,y;    printf("二分法:\n");    printf("误差          收敛度\n");    x1=a;    x2=b;    x3 = (x1+x2)/2;    while(fabs(x2-x1)>=zero)    {        x[i++]=x3;        //printf("%lf %lf %lf\n",x1,x2,x3);        if(fr(x1)*fr(x3)>0 )x1=x3;        else if(fr(x1)*fr(x3)<0)x2=x3;        else        {            printf("方程的根为%.7lf\n",x3);            break;        }        x3 = (x1+x2)/2;    }    for(j=0;j
        
         1)  短路特性    }    printf("\n经过%d运算,方程的根为%.7lf\n",i,x3);    return 0;}int diedai(double a){    int i=0,j;    double x1,x2;    printf("\n迭代法:\n");    printf("误差          收敛度\n");    //x1=a+1;                     //x1的初始化    x2=a;    do    {        x[i++]=x2;        x1 = x2;        //printf("%.7lf\n",x1);        x2 = fd(x1);    }while(fabs(x1-x2)>=zero);    for(j=0;j
         
          =zero);    for(j=0;j
         
        
       
      

#include
      
       #include
       
        #define zero 1e-5                   //允许的误差范围#define chuzhi 0.5double x[10000];                  //保存每一次的计算值double fr(double x){    return x-exp(-x);             //测试函数f(x)}double fd(double x){    return exp(-x);               //迭代法迭代式}double fn(double x)              //牛顿法迭代公式{    return x-(x-exp(-x))/(1+x);}int erfen(double a,double b){    int i=0,j;    double x1,x2,x3,y;    x1=a;    x2=b;    x3 = (x1+x2)/2;    while(fabs(x2-x1)>=chuzhi)    {        x[i++]=x3;        //printf("%lf %lf %lf\n",x1,x2,x3);        if(fr(x1)*fr(x3)>0 )x1=x3;        else if(fr(x1)*fr(x3)<0)x2=x3;        else        {            printf("方程的根为%.7lf\n",x3);            break;        }        x3 = (x1+x2)/2;    }    for(j=0;j
        
         =zero);    for(j=0;j