Question

求C语言课程设计：用高斯列主元消元法解线性方程组

运用所学的C语言，编制和调试程序，具有以下功能：用高斯列主元消元法解下列线性方程组：
2X1＋2X2＋5X3＝5
3X1＋4X2＋7X3＝6
X1＋3X2＋3X3＝5
并写下每个部骤的注释，和对编程的总体阐述！！！！急！！！！！
要每个部骤的注释，和对编程的总体阐述，因为还要答辩呢！！！！望高手指点！！！！

Answer 1

这里向你推荐一下克鲁特算法（其实就是对高斯列主元消元法进行优化，使之更适合于计算机编程），首先将矩阵A进行LU分解（将系数矩阵分解成一个上三角矩阵和一个下三角矩阵），分解的过程中用到了隐式的主元寻找法，同时利用克鲁特算法可以将两个n*n矩阵压缩到一个n*n矩阵中，大大节省了存储空间提高了计算速度。
方程可化为L*U*x=B，令U*x=y --->L*y=B
然后利用回代先求y，再利用y求x
因为该方法在求解过程中不涉及增广矩阵所以矩阵B几乎不参与什么运算，所以它的计算速度应该能够达到高斯列主元消元法的三倍，但原理与其基本一致。
而且我在程序中使用了动态数组方便你今后进行扩展。

以下程序按照《矩阵论第二版》和《C语言数值计算法方法大全》编写，LU分解部分程序主要参考了《C语言数值计算法方法大全》第二章的程序
如果你需要详细的理论讲解我可以将这两本书和源程序发给你，上面的论述相当详细足够你答辩用的了，我的邮箱hu_hu605@163.com

计算结果：
A矩阵:
2 2 5
3 4 7
1 3 3
B矩阵:
5
6
5
解矩阵:
x 1=-7
x 2=0.333333
x 3=3.66667

Press any key to continue

#include <cmath>
#include <iostream>
#include <iomanip>
#include <cstdlib>
#include <functional>
#include <vector>
#include <algorithm>
using namespace std;
#define TINY 1.0e-20 //A small number.
#define N 3
void ludcmp(vector<vector<float> > &a, int n, vector<int> &indx, float &d);//对矩阵进行LU分解
void lubksb(vector<vector<float> > &a, int n, vector<int> &indx, vector<float> &b);//对矩阵进行前向和后向回代
void root(vector<vector<float> > &x,vector<float> &col);//解方程结果保存在y中
void iniv(vector<vector<float> > &x,vector<float> line,int n);//对二维动态数组进行初始化
void main()
{
int i,j,n=N;//输入矩阵的维数
float A[N][N]={{2,2,5},{3,4,7},{1,3,3}};//左边A矩阵
float B[N]={5,6,5};//右边B矩阵
vector<vector<float> > x;//建立动态二维数组存放A，保证你的程序进行扩展时只改A，B，N
vector<float> line;
vector<float> y(n);//建立动态数组存放B
iniv(x,line,n);
y.clear();
for(i=0;i<n;i++)//将A赋给x，B赋给y
{
y.push_back(B[i]);
for(j=0;j<n;j++)
{
x[i].push_back(A[i][j]);
}
}
cout<<"A矩阵:"<<endl;
for(i=0;i<n;i++)
{
for(j=0;j<n;j++)
{
cout<<setw(2)<<setiosflags(ios::left)<<setw(2)<<x[i][j]<<" ";
}
cout<<endl;
}
cout<<"B矩阵:"<<endl;
for(i=0;i<n;i++)
{
cout<<setw(2)<<setiosflags(ios::left)<<setw(2)<<y[i]<<endl;
}
root(x,y);//求根
cout<<"解矩阵:"<<endl;
for(i=0;i<n;i++)
{
cout<<setw(2)<<setiosflags(ios::left)<<"x"<<i+1<<"="<<y[i]<<endl;
}
cout<<endl;
}
void root(vector<vector<float> > &x,vector<float> &col)
{
int n=x.size(),i=0,j=0;
vector<int> index(n);//用于记录寻找主元素过程中对矩阵的初等变换
index.clear();
float m=1.0;//记录变换方式，此程序中无用
ludcmp(x,n,index,m);//进行LU分解
lubksb(x,n,index,col);//根据分解结果进行回带
}
//以下程序按照《矩阵论第二版》和《C语言数值计算法方法大全》编写，LU分解部分程序主要参考了《C语言数值计算法方法大全》第二章的程序
//如果你需要详细的理论讲解我可以将这两本书和源程序发给你，我的邮箱hu_hu605@163.com
void ludcmp(vector<vector<float> > &a, int n, vector<int> &indx, float &d)
{

int i,imax,j,k;
float big=0,dum=0,sum=0,temp=0;
vector<float> vv(n);
vv.clear();
d=1.0;
for (i=0;i<n;i++)
{
big=0.0;
for (j=0;j<n;j++)
if ((temp=fabs(a[i][j])) > big)
big=temp;
vv[i]=1.0/big;
}
for (j=0;j<n;j++)
{
for (i=0;i<j;i++)
{
sum=a[i][j];
for (k=0;k<i;k++)
sum -= a[i][k]*a[k][j];
a[i][j]=sum;
}
big=0.0;
for (i=j;i<n;i++)
{
sum=a[i][j];
for (k=0;k<j;k++)
sum -= a[i][k]*a[k][j];
a[i][j]=sum;
if ( (dum=vv[i]*fabs(sum)) >= big)
{
big=dum;
imax=i;
}
}
if (j != imax)
{
for (k=0;k<n;k++)
{
dum=a[imax][k];
a[imax][k]=a[j][k];
a[j][k]=dum;
}
d = -(d);
vv[imax]=vv[j];
}
indx[j]=imax;
if (a[j][j] == 0.0)
a[j][j]=TINY;
if (j != n)
{
dum=1.0/(a[j][j]);
for (i=j+1;i<n;i++)
a[i][j] *= dum;
}
}
}
void lubksb(vector<vector<float> > &a, int n, vector<int> &indx, vector<float> &b)
{
int i,ii=0,ip,j;
float sum;
for(i=0;i<n;i++)//按LU分解时寻找主元所进行的初等变换进行反边变换。
{
ip=indx[i];
sum=b[ip];
b[ip]=b[i];
b[i]=sum;
}
sum=0;
for (i=1;i<n;i++)
{
sum=0;
for(j=0;j<i;j++)
{
sum+=a[i][j]*b[j];
}
b[i]=b[i]-sum;
}
b[n-1]=b[n-1]/a[n-1][n-1];
for (i=n-2;i>=0;i--)
{
sum=0;
for(j=i+1;j<n;j++)
{
sum+=a[i][j]*b[j];
}
b[i]=(b[i]-sum)/a[i][i];
}
}
void iniv(vector<vector<float> > &x,vector<float> line,int n)
{
int i,j;
for(i=0;i<n;i++)
{
x.push_back(line);
for(j=0;j<n;j++)
{
x[i].clear();
}
}
}

Answer 2

我在上面那位师兄提供的代码的基础上，稍稍做了一点修改，运行结果不仅可以显示X的值，还能显示出消元结束后系数矩阵转变成的上三角阵。代码如下：
#include<stdio.h>
#include<stdlib.h>
#include<conio.h>
#include<math.h>
#define n 4 /*方程组的阶数*/
#define precision 1e-16
static double aa[n][n+1]={{7.2,2.3,-4.4,0.5,15.1},{1.3,6.3,-3.5,\
2.8,1.8},{5.6,0.9,8.1,-1.3,16.6},{1.5,0.4,3.7,5.9,36.9}};
/*增广矩阵的原始数据*/
void main()
{int i,j,det;double a[n+1][n+2],x[n+1];
int GaussElimination_ColumnSelect();
clrscr();
for(i=1;i<=n;i++) for(j=1;j<=n+1;j++)
/*用a[1][1]~a[n][n+1]存放增广矩阵*/
a[i][j]=aa[i-1][j-1];
det=GaussElimination_ColumnSelect(a,x);
/*调用求解方程的函数，获取返回标志值*/
if(det!=0)
for(i=1;i<=n;i++)
printf("\nx[%d]=%f\n",i,x[i]);printf("\n");
getch();
}
int GaussElimination_ColumnSelect(double a[][n+2],double x[n+1])
/*用列主元高斯消去法求解线性方程组的函数*/
{int i,j,k,r;
double c;
for(k=1;k<=n-1;k++) /*消元过程*/
{r=k;
for(i=k;i<=n;i++) /*选列主元*/
if(fabs(a[i][k])>fabs(a[r][k]))r=i;
if(fabs(a[r][k])<precision)
{printf("\n det A=0.Elimination faild!");exit(0);}
if(r!=k)
{for(j=k;j<=n+1;j++) /*交换k,r两行*/
{c=a[k][j];a[k][j]=a[r][j];a[r][j]=c;}
}
for(i=k+1;i<=n;i++) /*进行消元计算*/
{c=a[i][k]/a[k][k];
for(j=k;j<=n+1;j++)
a[i][j]=a[i][j]-c*a[k][j];
}
}

printf("After modified,the matrix is::\n");
for(i=1;i<=n;i++){
for(j=1;j<=n;j++){
printf("%f\t",a[i][j]);
if(j%n==0) printf("\n");
}
}

if(fabs(a[n][n])<precision)
{printf("\n det A=0.Algorithm faild!");exit(0);}
for(k=n;k>=1;k--) /*回代过程*/
{x[k]=a[k][n+1];
for(j=k+1;j<=n;j++)
x[k]=x[k]-a[k][j]*x[j];
x[k]=x[k]/a[k][k];
}
return(1);
}

Answer 3

#include<stdio.h>
#include<stdlib.h>
#include<conio.h>
#include<math.h>
#define n 3 /*方程组的阶数*/
#define precision 1e-16
static double aa[n][n+1]={{2,2,5,5},{3,4,7,6},{1,3,3,5}};
/*增广矩阵的原始数据*/
void main()
{int i,j,det;double a[n+1][n+2],x[n+1];
int GaussElimination_ColumnSelect();

for(i=1;i<=n;i++) for(j=1;j<=n+1;j++)
/*用a[1][1]~a[n][n+1]存放增广矩阵*/
a[i][j]=aa[i-1][j-1];
det=GaussElimination_ColumnSelect(a,x);
/*调用求解方程的函数，获取返回标志值*/
if(det!=0)
for(i=1;i<=n;i++)
printf("\nx[%d]=%f\n",i,x[i]);printf("\n");
getch();
}
int GaussElimination_ColumnSelect(double a[][n+2],double x[n+1])
/*用列主元高斯消去法求解线性方程组的函数*/
{int i,j,k,r;double c;
for(k=1;k<=n-1;k++) /*消元过程*/
{r=k;
for(i=k;i<=n;i++) /*选列主元*/
if(fabs(a[i][k])>fabs(a[r][k]))r=i;
if(fabs(a[r][k])<precision)
{printf("\n det A=0.Elimination faild!");exit(0);}
if(r!=k)
{for(j=k;j<=n+1;j++) /*交换k,r两行*/
{c=a[k][j];a[k][j]=a[r][j];a[r][j]=c;}
}
for(i=k+1;i<=n;i++) /*进行消元计算*/
{c=a[i][k]/a[k][k];
for(j=k+1;j<=n+1;j++)
a[i][j]=a[i][j]-c*a[k][j];
}
}
if(fabs(a[n][n])<precision)
{printf("\n det A=0.Algorithm faild!");exit(0);}
for(k=n;k>=1;k--) /*回代过程*/
{x[k]=a[k][n+1];
for(j=k+1;j<=n;j++)
x[k]=x[k]-a[k][j]*x[j];
x[k]=x[k]/a[k][k];
}
return(1);
}

具体思想计算方法上有，在这也说不清。

Answer 4